#if defined(ROW_LAND)
#define SEA_P .true.
#define SEA_U .true.
#define SEA_V .true.
#elif defined(ROW_ALLSEA)
#define SEA_P allip(j).or.ip(i,j).ne.0
#define SEA_U alliu(j).or.iu(i,j).ne.0
#define SEA_V alliv(j).or.iv(i,j).ne.0
#else
#define SEA_P ip(i,j).ne.0
#define SEA_U iu(i,j).ne.0
#define SEA_V iv(i,j).ne.0
#endif
      module mod_tsadvc
      use mod_xc    ! HYCOM communication interface
      use mod_pipe  ! HYCOM debugging interface
      implicit none
!
! --- module for tsadvc and related routines
!
      private !! default is private
      public  :: tsadvc
!
      integer, parameter, dimension (0:4) :: &
        mbdy_advtyp = (/ 2,     & !PCM
                         5,     & !MPDATA
                         5,     & !FCT2
                         0,     & !N/A
                         5 /)     !FCT4
!
      logical, parameter :: lpipe_advem=.false.  !extra checking (when pipe on)
      logical, parameter :: lconserve  =.false.  !explicitly conserve the field
!
      logical, save ::  ldebug_tsdif   !switch to debug diffusion - usually .false.
      logical, save ::  ldebug_advem   !switch to debug advection - usually .false.
      integer, save ::  itests,jtests  !local copy of itest,jtest
!
#if defined(RELO)
      real, save, allocatable, dimension(:,:) :: &
#else
      real, save, dimension(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: &
#endif
         fmx,fmn       & ! local max,min
        ,flx,fly       & ! fluxes
        ,fldlo         & ! lo order solution
        ,fmxlo,fmnlo   & ! local min
        ,fax,fay       & ! fluxes
        ,rp,rm         & ! FCT/MPDATA terms
        ,flxdiv        & ! flux divergence
        ,tx1,ty1       & ! MPDATA terms
        ,fldao,fldan     ! total field quantity (old/center, new)

#if defined(RELO)
      real, save, allocatable, dimension(:,:) :: &
#else
      real,    save, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy):: &
#endif
          uloc,vloc,hloc,dtloc,ucumdt,vcumdt,flxcum,flycum

#if defined(RELO)
      logical, save, allocatable, dimension(:,:) :: &
#else
      logical, save, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy):: &
#endif
          lcalc

      contains

      subroutine advem(advtyp,fld,fldc,u,v,fco,fcn,posdef, &
                       scal,scali,dt2,btrmas)
      implicit none
!
      logical, intent(in)    :: btrmas
      integer, intent(in)    :: advtyp
      real,    intent(in)    :: posdef,dt2
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(inout) :: fld
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(in)    :: fldc,u,v,fco,fcn,scal,scali
!
! --- wrapper for advection schemes
!
! --- a recent text on advection schemes is:
! --- D.R. Durran (1999): Numerical Methods for wave equations in
! --- geophysical fluid dynamics, Springer.
!
! --- advtyp= 0 for 1st order PCM    (Donor Cell)
! --- advtyp= 1 for 2nd order MPDATA (old to new, as in 2.1.03)
! --- advtyp= 2 for 2nd order FCT    (Leapfrog time step)
! --- advtyp= 4 for 4th order FCT    (Leapfrog time step)
!
! --- time steps are "old", "center" and "new".
!
! --- fld    - scalar field, at old time step on input but new on output
! --- fldc   - scalar field, at center time step
! --- u,v    - mass fluxes satisfying continuity equation (old to new)
! --- fco    - thickness of the layer at old time step
! --- fcn    - thickness of the layer at new time step
! --- posdef - offset for MPDATA to make the field positive
! --- scal   - spatial increments (squared)
! --- scali  - inverse of scal
! --- dt2    - temporal increment (from old to new, i.e. two time steps)
!
!  on return, fld's valid halo will be 0 wide.
!
      real    offset
      real*8  sumold,sumnew,sumcor
      integer i,j
!
#if defined(RELO)
      if     (.not.allocated(fmx)) then
        allocate( &
                 fmx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 fmn(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 flx(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 fly(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               fldlo(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               fmxlo(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               fmnlo(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 fax(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 fay(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                  rp(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                  rm(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
              flxdiv(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 tx1(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 ty1(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               fldao(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               fldan(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) )
        call mem_stat_add( 16*(idm+2*nbdy)*(jdm+2*nbdy) )
                 fmx = r_init
                 fmn = r_init
                 flx = r_init
                 fly = r_init
               fldlo = r_init
               fmxlo = r_init
               fmnlo = r_init
                 fax = r_init
                 fay = r_init
                  rp = r_init
                  rm = r_init
              flxdiv = r_init
                 tx1 = r_init
                 ty1 = r_init
               fldao = r_init
               fldan = r_init
      endif !.not.allocated
#endif
!
      if     (advtyp.eq.0) then
        call advem_pcm(   fld,     u,v,fco,fcn,       scal,scali,dt2)
      elseif (advtyp.eq.1) then
        call advem_mpdata(fld,     u,v,fco,fcn,posdef,scal,scali,dt2)
      elseif (advtyp.eq.2 .and. btrmas) then
        call advem_fct2c( fld,fldc,u,v,fco,fcn,       scal,scali,dt2)
      elseif (advtyp.eq.2) then !.not.btrmas
        call advem_fct2(  fld,fldc,u,v,fco,fcn,       scal,scali,dt2)
      elseif (advtyp.eq.4) then
        call advem_fct4(  fld,fldc,u,v,fco,fcn,       scal,scali,dt2)
      else
        if     (mnproc.eq.1) then
          write(lp,'(/ a,i4 /)') &
            'error: advem called with advtyp =',advtyp
        endif
        call xcstop('advem')
               stop 'advem'
      endif
!
      if     (lconserve) then !usually .false.
!
! ---   explicit conservation of tracer (should not be needed).
!
        call xcsum(sumold, fldao,ipa)
        call xcsum(sumnew, fldan,ipa)
!
        if     (sumnew.ne.0.0) then
          offset = (sumold-sumnew)/sumnew
        else
          offset = 0.0
        endif
!
!$OMP   PARALLEL DO PRIVATE(j,i,offset) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1,jj
          do i=1,ii
            if (SEA_P) then
              fld(i,j)=fld(i,j)*(1.0+offset)
!
!diag         fldan(i,j) = fld(i,j)*fcn(i,j)*scal(i,j)
            endif !ip
          enddo !i
        enddo !j
!$OMP   END PARALLEL DO
        if     (lpipe .and. lpipe_advem) then
! ---     compare two model runs.
          call pipe_compare_sym1(fld,    ip,'ad:oset:fld ')
        endif
!
!diag   call xcsum(sumcor, fldan,ipa)
!diag   if     (mnproc.eq.1) then
!diag     write(lp,'(a,1p4e16.8)') &
!diag       'advem: ',sumold,sumnew,sumcor,offset
!diag   endif
      endif !lconserve
      return
      end subroutine advem

      subroutine advem_mpdata(fld,u,v,fco,fcn,posdef,scal,scali,dt2)
      implicit none
!
      real,    intent(in)    :: posdef,dt2
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(inout) :: fld
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(in)    :: u,v,scal,scali,fco,fcn
!
! LeapFrog 2nd order MPDATA.
! combined monotone scheme, for details see section 3.3 (eqs. 34 to 37)
! in smolarkiewicz and clark, 1986, j.comput.phys.,67,no 2, p. 396-438
! and smolarkiewicz and grabowski, 1989, j.comput.phys. and recently
! P.K. Smolarkiewicz and L.J. Margolin (1998): MPDATA: A finite
! difference solver for geophysical flows, J.Comput.Phys. 140 459-480.
!
! time steps are "old", "center" and "new".
!
!  fld    - scalar field, must be >0, old input but new output
!  u,v    - mass fluxes satisfying continuity equation (old to new)
!  fco    - thickness of the layer at old time step
!  fcn    - thickness of the layer at new time step
!  posdef - offset to make the field positive
!  scal   - spatial increments (squared)
!  scali  - inverse of scal
!  dt2    - temporal increment (from old to new)
!
!  on return, fld's valid halo will be 0 wide.
!
      real,    parameter :: onemu=9806.e-12 !very small layer thickness
!
      real    fcn2,fco2,flxdn,flxdp,flydn,flydp,q
      integer i,j,l,ia,ib,ja,jb,mbdy_a,margin
!
      mbdy_a = mbdy_advtyp(1)  ! = 5
!
! --- compute low-order and part of antidiffusive fluxes
!
! --- rhs: u, v, fld+
! --- lhs: flx, fly, fmx, fmn
!
      margin = mbdy_a - 1
!
!$OMP PARALLEL DO PRIVATE(j,i,q,jb,ja,ib,ia) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            tx1(i,j)=.5*abs(u(i,j))*(fld(i,j)-fld(i-1,j))
            if (u(i,j).ge.0.0) then
              q=fld(i-1,j)
            else
              q=fld(i  ,j)
            endif
            flx(i,j)=u(i,j)*(q+posdef)
          endif !iu
          if (SEA_V) then
            ty1(i,j)=.5*abs(v(i,j))*(fld(i,j)-fld(i,j-1))
            if (v(i,j).ge.0.0) then
              q=fld(i,j-1)
            else
              q=fld(i,j  )
            endif
            fly(i,j)=v(i,j)*(q+posdef)
          endif !iv
          if (SEA_P) then
            ia=ipim1(i,j) !i-1 if sea and i if land
            ib=ipip1(i,j) !i+1 if sea and i if land
            ja=ipjm1(i,j) !j-1 if sea and j if land
            jb=ipjp1(i,j) !j+1 if sea and j if land
            fmx(i,j)=max(fld(i,j),fld(ia,j),fld(ib,j), &
                                  fld(i,ja),fld(i,jb))+posdef
            fmn(i,j)=min(fld(i,j),fld(ia,j),fld(ib,j), &
                                  fld(i,ja),fld(i,jb))+posdef
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(tx1, iu,'ad:11:tx1   ', &
                               ty1, iv,'ad:11:ty1   ')
        call pipe_compare_sym2(flx, iu,'ad:11:flx   ', &
                               fly, iv,'ad:11:fly   ')
        call pipe_compare_sym1(fmx, ip,'ad:11:fmx   ')
        call pipe_compare_sym1(fmn, ip,'ad:11:fmn   ')
      endif
!
! --- rhs: u, v, fld+
! --- lhs: flx, fly, fmx, fmn
!
      margin = mbdy_a - 1
!
      do j=1-margin,jj+margin
        do l=1,isp(j) !ok
          if     (ifp(j,l).ge. 1-margin) then
            flx(ifp(j,l)  ,j)=0.0
          endif
          if     (ilp(j,l).lt.ii+margin) then
            flx(ilp(j,l)+1,j)=0.0
          endif
        enddo !l
      enddo !j
!
      do i=1-margin,ii+margin
        do l=1,jsp(i) !ok
          if     (jfp(i,l).ge. 1-margin) then
            fly(i,jfp(i,l)  )=0.0
          endif
          if     (jlp(i,l).lt.jj+margin) then
            fly(i,jlp(i,l)+1)=0.0
          endif
        enddo !l
      enddo !i
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(flx, iu,'ad:22:flx   ', &
                               fly, iv,'ad:33:fly   ')
      endif
!diag if     (itests.gt.0 .and. jtests.gt.0) then
!diag   i=itests
!diag   j=jtests
!diag   write (lp,'(a,2i5,f22.3/1pe39.2/0pf21.3,1pe9.2,0pf9.3, &
!diag              1pe9.2,0pf9.3/1pe39.2/0pf39.3)') &
!diag     'advem (1)',i+i0,j+j0, &
!diag     fld(i-1,j),u(i,j),fld(i,j-1),v(i,j), &
!diag     fld(i,j),v(i,j+1),fld(i,j+1),u(i+1,j),fld(i+1,j)
!diag endif
!
! --- rhs: flx+, fly+, fco, fmn, fmx, fcn
! --- lhs: fld
!
      margin = mbdy_a - 2
!
!$OMP PARALLEL DO PRIVATE(j,i,q) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
!...........lo order Donor Cell step
            flxdiv(i,j)=((flx(i+1,j)-flx(i,j))+ &
                         (fly(i,j+1)-fly(i,j)) )*dt2*scali(i,j)
            q=(fld(i,j)+posdef)*(fco(i,j)+onemu)-flxdiv(i,j)
            !max,min should only be active for very thin layers
            fldlo(i,j)=max( fmn(i,j), min( fmx(i,j), &
                                           q/(fcn(i,j)+onemu) ))
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym1(fldlo,  ip,'ad:610:fldlo')
        call pipe_compare_sym1(flxdiv, ip,'ad:610:flxdv')
      endif
!
! --- finish computation of antidiffusive fluxes
!
! --- rhs: tx1, u, ty1, v, flxdiv+, fco+, fcn+
! --- lhs: flx, fly
!
      margin = mbdy_a - 2
!
!$OMP PARALLEL DO PRIVATE(j,i,fcn2,fco2) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        fco2=0.0
        fcn2=0.0
        do i=1-margin,ii+margin
          if (SEA_U) then
            fco2=fco(i,j)+fco(i-1,j)  ! inforce order on flx calc
            fcn2=fcn(i,j)+fcn(i-1,j)  ! inforce order on flx calc
            flx(i,j)=tx1(i,j)-u(i,j)*(flxdiv(i,j)+flxdiv(i-1,j)) &
               /((fco2+fcn2)+onemu)
          endif !iu
          if (SEA_V) then
            fco2=fco(i,j)+fco(i,j-1)  ! inforce order on fly calc
            fcn2=fcn(i,j)+fcn(i,j-1)  ! inforce order on fly calc
            fly(i,j)=ty1(i,j)-v(i,j)*(flxdiv(i,j)+flxdiv(i,j-1)) &
               /((fco2+fcn2)+onemu)
          endif !iv
        enddo !i
        if (fco2*fcn2.eq.1.e30) flx(1-nbdy,j)=0.0  ! prevent removal of fc*2
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(flx,  iu,'ad: 8:flx   ', &
                               fly,  iv,'ad: 8:fly   ')
      endif
!
! --- limit antidiffusive fluxes
! --- rp and rm used to be called flp and fln
!
! --- rhs: fmx, fmn, fldlo, fcn, flx+, fly+
! --- lhs: rp, rm
!
      margin = mbdy_a - 3
!
!$OMP PARALLEL DO PRIVATE(j,i,flxdn,flxdp,flydn,flydp) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
            flxdp=min(0.0,flx(i+1,j))-max(0.0,flx(i,j))
            flxdn=max(0.0,flx(i+1,j))-min(0.0,flx(i,j))
            flydp=min(0.0,fly(i,j+1))-max(0.0,fly(i,j))
            flydn=max(0.0,fly(i,j+1))-min(0.0,fly(i,j))
            rp(i,j)=(fmx(i,j)-fldlo(i,j))*(fcn(i,j)*scal(i,j))/ &
             ((onemu-(flxdp+flydp))*dt2)
            rm(i,j)=(fldlo(i,j)-fmn(i,j))*(fcn(i,j)*scal(i,j))/ &
             ((onemu+(flxdn+flydn))*dt2)
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym1(rp, ip,'ad:16:flp   ')
        call pipe_compare_sym1(rm, ip,'ad:16:fln   ')
      endif
!
! --- rhs: flx, fly, rp+, rm+
! --- lhs: flx, fly
!
      margin = mbdy_a - 4
!
!$OMP PARALLEL DO PRIVATE(j,i) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            flx(i,j)=max(0.0,flx(i,j))*min(1.0,rp(i,j),rm(i-1,j)) &
                    +min(0.0,flx(i,j))*min(1.0,rp(i-1,j),rm(i,j))
          endif !iu
          if (SEA_V) then
            fly(i,j)=max(0.0,fly(i,j))*min(1.0,rp(i,j),rm(i,j-1)) &
                    +min(0.0,fly(i,j))*min(1.0,rp(i,j-1),rm(i,j))
          endif !iv
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(flx, iu,'ad:18:flx   ', &
                               fly, iv,'ad:18:fly   ')
      endif
!
!diag i=itests
!diag j=jtests
!diag write (lp,'(''advem (2)'',2i5,f22.3/1pe39.2/0pf21.3,1pe9.2,0pf9.3, &
!diag 1pe9.2,0pf9.3/1pe39.2/0pf39.3)') i,j,fldlo(i-1,j),u(i,j),fldlo(i,j-1), &
!diag v(i,j),fldlo(i,j),v(i,j+1),fldlo(i,j+1),u(i+1,j),fldlo(i+1,j)
!
! --- rhs: flx+, fly+, fldlo, fcn, fmx
! --- lhs: flxdiv, fld
!
      margin = mbdy_a - 5
!
!$OMP PARALLEL DO PRIVATE(j,i) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
            fldao(i,j) = fld(i,j)*fco(i,j)*scal(i,j)
!
!...........apply antidiffusive flux correction
            flxdiv(i,j)=((flx(i+1,j)-flx(i,j))+ &
                         (fly(i,j+1)-fly(i,j)) )*dt2*scali(i,j)
            !max,min should only be active for very thin layers
            fld(i,j)=max( fmn(i,j), min( fmx(i,j), &
                          fldlo(i,j)-flxdiv(i,j)/(fcn(i,j)+onemu) ))
            fld(i,j)=fld(i,j)-posdef
!
            fldan(i,j) = fld(i,j)*fcn(i,j)*scal(i,j)
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym1(flxdiv, ip,'ad:620:flxdv')
        call pipe_compare_sym1(fld,    ip,'ad:1620:fld ')
      endif
      return
      end subroutine advem_mpdata

      subroutine advem_pcm(fld,u,v,fco,fcn,scal,scali,dt2)
      implicit none
!
      real,    intent(in)    :: dt2
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(inout) :: fld
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(in)    :: u,v,scal,scali,fco,fcn
!
! Piecewise Constant Method (Donor Cell, Upwind)
! Over two time steps (may require half the normal time step for stability).
!
! time steps are "old", "center" and "new".
!
!  fld   - scalar field, need not be >0, old input but new output
!  u,v   - mass fluxes satisfying continuity equation (old to new)
!  fco   - thickness of the layer at old    time step
!  fcn   - thickness of the layer at new    time step
!  scal  - spatial increments (squared)
!  scali - inverse of scal
!  dt2   - temporal increment (from old to new)
!
!  on return, fld's valid halo will be 0 wide.
!
      real,    parameter :: onemu=9806.e-12 !very small layer thickness
!
      real    q
      integer i,j,l,ia,ib,ja,jb,mbdy_a,margin
!
! --- rhs: u, v, fld+
! --- lhs: flx, fly
!
      mbdy_a = mbdy_advtyp(0)  ! = 2
!
      margin = mbdy_a - 1
!
!$OMP PARALLEL DO PRIVATE(j,i,q) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            if (u(i,j).ge.0.0) then
              q=fld(i-1,j)
            else
              q=fld(i  ,j)
            endif
            flx(i,j)=u(i,j)*q
          endif !iu
          if (SEA_V) then
            if (v(i,j).ge.0.0) then
              q=fld(i,j-1)
            else
              q=fld(i,j  )
            endif
            fly(i,j)=v(i,j)*q
          endif !iv
          if (SEA_P) then
            ia=ipim1(i,j) !i-1 if sea and i if land
            ib=ipip1(i,j) !i+1 if sea and i if land
            ja=ipjm1(i,j) !j-1 if sea and j if land
            jb=ipjp1(i,j) !j+1 if sea and j if land
            fmx(i,j)=max(fld(i,j),fld(ia,j),fld(ib,j), &
                                  fld(i,ja),fld(i,jb))
            fmn(i,j)=min(fld(i,j),fld(ia,j),fld(ib,j), &
                                  fld(i,ja),fld(i,jb))
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(flx, iu,'ad:11:flx   ', &
                               fly, iv,'ad:11:fly   ')
      endif
!
      do j=1-margin,jj+margin
        do l=1,isp(j) !ok
          if     (ifp(j,l).ge. 1-margin) then
            flx(ifp(j,l)  ,j)=0.0
          endif
          if     (ilp(j,l).lt.ii+margin) then
            flx(ilp(j,l)+1,j)=0.0
          endif
        enddo !l
      enddo !j
!
      do i=1-margin,ii+margin
        do l=1,jsp(i) !ok
          if     (jfp(i,l).ge. 1-margin) then
            fly(i,jfp(i,l)  )=0.0
          endif
          if     (jlp(i,l).lt.jj+margin) then
            fly(i,jlp(i,l)+1)=0.0
          endif
        enddo !l
      enddo !i
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(flx, iu,'ad:22:flx   ', &
                               fly, iv,'ad:33:fly   ')
      endif
!diag if     (itests.gt.0 .and. jtests.gt.0) then
!diag   i=itests
!diag   j=jtests
!diag   write (lp,'(a,2i5,f22.3/1pe39.2/0pf21.3,1pe9.2,0pf9.3, &
!diag              1pe9.2,0pf9.3/1pe39.2/0pf39.3)') &
!diag     'advem (1)',i+i0,j+j0, &
!diag     fld(i-1,j),u(i,j),fld(i,j-1),v(i,j), &
!diag     fld(i,j),v(i,j+1),fld(i,j+1),u(i+1,j),fld(i+1,j)
!diag endif
!
! --- rhs: flx+, fly+, fld, fco, fcn
! --- lhs: flxdiv, fld
!
      margin = mbdy_a - 2
!
!$OMP PARALLEL DO PRIVATE(j,i,q) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
            fldao(i,j) = fld(i,j)*fco(i,j)*scal(i,j)
!
            flxdiv(i,j)=((flx(i+1,j)-flx(i,j))+ &
                         (fly(i,j+1)-fly(i,j)) )*dt2*scali(i,j)
            q=fld(i,j)*(fco(i,j)+onemu)-flxdiv(i,j)
            !max,min should only be active for very thin layers
            fld(i,j)=max( fmn(i,j), min( fmx(i,j), &
                                         q/(fcn(i,j)+onemu) ))
!
            fldan(i,j) = fld( i,j)*fcn(i,j)*scal(i,j)
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym1(fld,    ip,'ad:610:fld  ')
        call pipe_compare_sym1(flxdiv, ip,'ad:610:flxdv')
      endif
!
!diag i=itests
!diag j=jtests
!diag write (lp,'(''advem (2)'',2i5,f22.3/1pe39.2/0pf21.3,1pe9.2,0pf9.3, &
!diag 1pe9.2,0pf9.3/1pe39.2/0pf39.3)') i,j,fld(i-1,j),u(i,j),fld(i,j-1), &
!diag v(i,j),fld(i,j),v(i,j+1),fld(i,j+1),u(i+1,j),fld(i+1,j)
      return
      end subroutine advem_pcm

      subroutine advem_fct2(fld,fldc,u,v,fco,fcn,scal,scali,dt2)
      implicit none
!
      real,    intent(in)    :: dt2
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(inout) :: fld
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(in)    :: fldc,u,v,scal,scali,fco,fcn
!
! Leapfrog 2nd order FCT
! S.T. Zalesak (1979): Fully multidimensional flux-corrected
! transport algorithms for fluids, J.Comput.Phys. 31 335-362.
!
! time steps are "old", "center" and "new".
!
!  fld   - scalar field, need not be >0, old input but new output
!  fldc  - scalar field at center time step
!  u,v   - mass fluxes satisfying continuity equation (old to new)
!  fco   - thickness of the layer at old    time step
!  fcn   - thickness of the layer at new    time step
!  scal  - spatial increments (squared)
!  scali - inverse of scal
!  dt2   - temporal increment (from old to new)
!
!  on return, fld's valid halo will be 0 wide.
!
      real,    parameter :: onemu=9806.e-12 !very small layer thickness
!
      real    flxdn,flxdp,flydn,flydp,q
      integer i,j,l,ia,ib,ja,jb,mbdy_a,margin
!
      real :: fhx, fhy, fqmax, fqmin, famax, famin
      real :: qdt2, qp, qm, fact

      mbdy_a = mbdy_advtyp(2)  ! = 5
!
! --- compute low-order and part of antidiffusive fluxes
!
! --- rhs: u, v, fld+
! --- lhs: flx, fly, fmx, fmn
!
      margin = mbdy_a - 1
!
!$OMP PARALLEL DO PRIVATE(j,i,q,jb,ja,ib,ia) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            if (u(i,j).ge.0.0) then
              q=fld(i-1,j)
            else
              q=fld(i  ,j)
            endif
            flx(i,j)=u(i,j)*q
          endif !iu
          if (SEA_V) then
            if (v(i,j).ge.0.0) then
              q=fld(i,j-1)
            else
              q=fld(i,j  )
            endif
            fly(i,j)=v(i,j)*q
          endif !iv
          if (SEA_P) then
            ia=ipim1(i,j) !i-1 if sea and i if land
            ib=ipip1(i,j) !i+1 if sea and i if land
            ja=ipjm1(i,j) !j-1 if sea and j if land
            jb=ipjp1(i,j) !j+1 if sea and j if land
            fmx(i,j)=max(fld(i,j),fld(ia,j),fld(ib,j), &
                                  fld(i,ja),fld(i,jb))
            fmn(i,j)=min(fld(i,j),fld(ia,j),fld(ib,j), &
                                  fld(i,ja),fld(i,jb))
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
!       call pipe_compare_sym2(tx1, iu,'ad:11:tx1   ',
!    &                         ty1, iv,'ad:11:ty1   ')
!       call pipe_compare_sym2(flx, iu,'ad:11:flx   ',
!    &                         fly, iv,'ad:11:fly   ')
      endif
!
      if     (ldebug_advem .and. itests.gt.0 .and. jtests.gt.0) then
        i=itests
        j=jtests
        write (lp,'(a,2i5,3f10.5)') &
          'advem: fld, rng  ',i+i0,j+j0, &
          fld(i,j),fmx(i,j),fmn(i,j)
      endif
!
      do j=1-margin,jj+margin
        do l=1,isp(j) !ok
          if     (ifp(j,l).ge. 1-margin) then
            flx(ifp(j,l)  ,j)=0.0
          endif
          if     (ilp(j,l).lt.ii+margin) then
            flx(ilp(j,l)+1,j)=0.0
          endif
        enddo !l
      enddo !j
!
      do i=1-margin,ii+margin
        do l=1,jsp(i) !ok
          if     (jfp(i,l).ge. 1-margin) then
            fly(i,jfp(i,l)  )=0.0
          endif
          if     (jlp(i,l).lt.jj+margin) then
            fly(i,jlp(i,l)+1)=0.0
          endif
        enddo !l
      enddo !i
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(flx, iu,'ad:22:flx   ', &
                               fly, iv,'ad:33:fly   ')
      endif
!
!diag if     (ldebug_advem .and. itests.gt.0 .and. jtests.gt.0) then
!diag   i=itests
!diag   j=jtests
!diag   write (lp,'(a,2i5,f22.3/1pe39.2/0pf21.3,1pe9.2,0pf9.3, &
!diag              1pe9.2,0pf9.3/1pe39.2/0pf39.3)') &
!diag     'advem (1)',i+i0,j+j0, &
!diag     fld(i-1,j),u(i,j),fld(i,j-1),v(i,j), &
!diag     fld(i,j),v(i,j+1),fld(i,j+1),u(i+1,j),fld(i+1,j)
!diag endif
!
! --- rhs: flx+, fly+, fld, fmx, fmn, fco, fcn
! --- lhs: fldlo, fmxlo, fmnlo, flxdiv
!
      margin = mbdy_a - 2
!
!$OMP PARALLEL DO PRIVATE(j,i,q) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
!...........lo order Donor Cell step
            flxdiv(i,j)=((flx(i+1,j)-flx(i,j))+ &
                         (fly(i,j+1)-fly(i,j)) )*dt2*scali(i,j)
            q=fld(i,j)*(fco(i,j)+onemu)-flxdiv(i,j)
            !max,min should only be active for very thin layers
            fldlo(i,j)=max( fmn(i,j), min( fmx(i,j), &
                                           q/(fcn(i,j)+onemu) ))
            fmxlo(i,j) = max(fld(i,j),fldc(i,j),fldlo(i,j))
            fmnlo(i,j) = min(fld(i,j),fldc(i,j),fldlo(i,j))
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym1(fldlo,  ip,'ad:610:fldlo')
        call pipe_compare_sym1(flxdiv, ip,'ad:610:flxdv')
      endif
!
      if     (ldebug_advem .and. itests.gt.0 .and. jtests.gt.0) then
        i=itests
        j=jtests
        write (lp,'(a,2i5,3f10.5)') &
          'advem: fldlo,rng ',i+i0,j+j0, &
          fldlo(i,j),fmnlo(i,j),fmxlo(i,j)
      endif
!
!.....Leapfrog step using high order scheme
!
! --- rhs: u, v, fld+, flx, fly
! --- lhs: fax, fay
!
      margin = mbdy_a - 2
!
!$OMP PARALLEL DO PRIVATE(j,i,q,jb,ja,ib,ia,fhx,fhy) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            fhx=u(i,j)*0.5*(fldc(i,j)+fldc(i-1,j))  ! 2nd order in space
            fax(i,j)= fhx-flx(i,j)                  ! anti-diffusion x-flux
          endif !iu
          if (SEA_V) then
            fhy=v(i,j)*0.5*(fldc(i,j)+fldc(i,j-1))  ! 2nd order in space
            fay(i,j)= fhy-fly(i,j)                  ! anti-diffusion y-flux
          endif !iv
        enddo !i
      enddo !j
!$OMP END PARALLEL DO

      do j=1-margin,jj+margin
        do l=1,isp(j) !ok
          if     (ifp(j,l).ge. 1-margin) then
            fax(ifp(j,l)  ,j)=0.0
          endif
          if     (ilp(j,l).lt.ii+margin) then
            fax(ilp(j,l)+1,j)=0.0
          endif
        enddo !l
      enddo !j
!
      do i=1-margin,ii+margin
        do l=1,jsp(i) !ok
          if     (jfp(i,l).ge. 1-margin) then
            fay(i,jfp(i,l)  )=0.0
          endif
          if     (jlp(i,l).lt.jj+margin) then
            fay(i,jlp(i,l)+1)=0.0
          endif
        enddo !l
      enddo !i
!========================================================
!
! --- finish computation of antidiffusive fluxes
!
! --- rhs: fmnlo+,fmxlo+,fax+,fay+,fldlo,fcn,scal
! --- lhs: rp, rm
!
      margin = mbdy_a - 3
!
      qdt2 = 1.0/dt2
!$OMP PARALLEL DO PRIVATE(j,i,ia,ib,ja,jb, &
!$OMP                     fqmax,fqmin,famax,famin,qp,qm) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
            ia=ipim1(i,j) !i-1 if sea and i if land
            ib=ipip1(i,j) !i+1 if sea and i if land
            ja=ipjm1(i,j) !j-1 if sea and j if land
            jb=ipjp1(i,j) !j+1 if sea and j if land
            fqmax = max(fmxlo(i,j),fmxlo(ia,j),fmxlo(ib,j), &
                                   fmxlo(i,ja),fmxlo(i,jb))
            fqmin = min(fmnlo(i,j),fmnlo(ia,j),fmnlo(ib,j), &
                                   fmnlo(i,ja),fmnlo(i,jb))
            famax = max(0.0,fax(i, j) ) - min(0.0,fax(ib,j) ) + &
                    max(0.0,fay(i, j) ) - min(0.0,fay(i, jb))
            famin = max(0.0,fax(ib,j) ) - min(0.0,fax(i, j) ) + &
                    max(0.0,fay(i, jb)) - min(0.0,fay(i, j) )
            if (famax > 0.0) then
              qp = (fqmax-fldlo(i,j)) *fcn(i,j)*scal(i,j)*qdt2
              if (qp < famax) then
                rp(i,j) = qp/famax  !<1.0
              else
                rp(i,j) = 1.0
              endif
            else
              rp(i,j) = 0.0
            endif
            if (famin > 0.0) then
              qm = (fldlo(i,j)-fqmin) *fcn(i,j)*scal(i,j)*qdt2
              if (qm < famin) then
                rm(i,j) = qm/famin  !<1.0
              else
                rm(i,j) = 1.0
              endif
            else
              rm(i,j) = 0.0
            endif
            fmx(i,j) = fqmax  !less restrictive maximum
            fmn(i,j) = fqmin  !less restrictive minimum
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
      if     (ldebug_advem .and. itests.gt.0 .and. jtests.gt.0) then
        i=itests
        j=jtests
        write (lp,'(a,2i5,3f10.5)') &
          'advem: fldc,rngfq',i+i0,j+j0, &
          fldc(i,j),fmn(i,j),fmx(i,j)
      endif
!
! --- rhs: rp+, rm+
! --- lhs: fax, fay
!
      margin = mbdy_a - 4
!
!$OMP PARALLEL DO PRIVATE(j,i,fact) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            if (fax(i,j) < 0.0) then
              fact = min(rp(i-1,j),rm(i,j))
            else
              fact = min(rp(i,j),rm(i-1,j))
            endif
            fax(i,j) = fact*fax(i,j)
          endif !iu
          if (SEA_V) then
            if (fay(i,j) < 0.0) then
              fact = min(rp(i,j-1),rm(i,j))
            else
              fact = min(rp(i,j),rm(i,j-1))
            endif
            fay(i,j) = fact*fay(i,j)
          endif !iv
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(flx, iu,'ad:18:flx   ', &
                               fly, iv,'ad:18:fly   ')
      endif
!
!diag i=itests
!diag j=jtests
!diag write (lp,'(''advem (2)'',2i5,f22.3/1pe39.2/0pf21.3,1pe9.2,0pf9.3, &
!diag 1pe9.2,0pf9.3/1pe39.2/0pf39.3)') i,j,fldlo(i-1,j),u(i,j),fldlo(i,j-1), &
!diag v(i,j),fldlo(i,j),v(i,j+1),fldlo(i,j+1),u(i+1,j),fldlo(i+1,j)
!
! --- rhs: fax+, fay+, fld, fcn, fmx, fmn
! --- lhs: fld
!
      margin = mbdy_a - 5
!
!$OMP PARALLEL DO PRIVATE(j,i) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
            fldao(i,j) = fld(i,j)*fco(i,j)*scal(i,j)
!
!...........apply antidiffusive flux correction
            flxdiv(i,j)=((fax(i+1,j)-fax(i,j))+ &
                         (fay(i,j+1)-fay(i,j)) )*dt2*scali(i,j)
            !max,min should only be active for very thin layers
            fld(i,j)=max( fmn(i,j), min( fmx(i,j), &
                          fldlo(i,j)-flxdiv(i,j)/(fcn(i,j)+onemu) ))
!
            fldan(i,j) = fld(i,j)*fcn(i,j)*scal(i,j)
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym1(flxdiv, ip,'ad:620:flxdv')
        call pipe_compare_sym1(fld,    ip,'ad:1620:fld ')
      endif
!
      if     (ldebug_advem .and. itests.gt.0 .and. jtests.gt.0) then
        i=itests
        j=jtests
        write (lp,'(a,2i5,f10.5)') &
          'advem: fld       ',i+i0,j+j0, &
          fld(i,j)
      endif
      return
      end subroutine advem_fct2

      subroutine advem_fct2c(fld,fldc,u,v,fco,fcn, scal,scali,dt2)
      implicit none
!
      real,    intent(in)    :: dt2
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(inout) :: fld
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(in)    :: fldc,u,v,scal,scali,fco,fcn
!
! Leapfrog 2nd order FCT
! S.T. Zalesak (1979): Fully multidimensional flux-corrected
! transport algorithms for fluids, J.Comput.Phys. 31 335-362.
!
! time steps are "old", "center" and "new".
!
!  fld   - scalar field, need not be >0, old input but new output
!  fldc  - scalar field at center time step
!  u,v   - mass fluxes satisfying continuity equation (old to new)
!  fco   - thickness of the layer at old    time step
!  fcn   - thickness of the layer at new    time step
!  scal  - spatial increments (squared)
!  scali - inverse of scal
!  dt2   - temporal increment (from old to new)
!
!  on return, fld's valid halo will be 0 wide.
!
!  Remy Baraille, SHOM, September 2008.
!
      real,    parameter :: epsil=1.e-10
!
      real    flxdn,flxdp,flydn,flydp,q
      integer i,j,l,ia,ib,ja,jb,mbdy_a
      integer iter,margin
!
      real :: fhx, fhy, fqmax, fqmin, famax, famin
      real :: qdt2, qp, qm, fact
!
#if defined(RELO)
      if     (.not.allocated(uloc)) then
       allocate( &
                 uloc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 vloc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 hloc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                dtloc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               ucumdt(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               vcumdt(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               flxcum(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               flycum(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                lcalc(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy))
        call mem_stat_add( 9*(idm+2*nbdy)*(jdm+2*nbdy) )
                 uloc = r_init
                 vloc = r_init
                 hloc = r_init
                dtloc = r_init
               ucumdt = r_init
               vcumdt = r_init
               flxcum = r_init
               flycum = r_init
               lcalc  = .false.
      endif !.not.allocated
#endif
!
      mbdy_a = mbdy_advtyp(2)  ! = 5
!
! --- compute low-order and part of antidiffusive fluxes
!
! --- rhs: u, v, fld+
! --- lhs: flx, fly, fmx, fmn
!
       lcalc(:,:) = .true.
      flxcum(:,:) = 0.0
      flycum(:,:) = 0.0
       dtloc(:,:) = 0.0
        hloc(:,:) = fco(:,:)
       fldlo(:,:) = fld(:,:)
!
! --- explicitly set to zero so that arrays are zero over land
      ucumdt(:,:) = 0.0
      vcumdt(:,:) = 0.0
        uloc(:,:) = 0.0
        vloc(:,:) = 0.0
         flx(:,:) = 0.0
         fly(:,:) = 0.0
!
      do iter=1,5
!
        margin = mbdy_a
!
!$OMP   PARALLEL DO PRIVATE(j,i,q) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-margin, jj+margin
          do i=1-margin,ii+margin
            if (SEA_P) then
              q = max(u(i+1,j),0.0)-min(u(i,j),0.0) &
                + max(v(i,j+1),0.0)-min(v(i,j),0.0)
              if (q.gt.0.0) then
                 dtloc(i,j) = min( dt2, hloc(i,j) / (q*scali(i,j)))
              else
                 dtloc(i,j) = dt2
              endif
            endif !ip
          enddo !i
        enddo !j
!
        margin = mbdy_a - 1
!
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-margin,jj+margin
          do i=1-margin,ii+margin
            if (SEA_U) then
              if (ucumdt(i,j).ne.dt2) then
                if (u(i,j).ge.0) then
                  uloc(i,j) = min(dt2-ucumdt(i,j),dtloc(i-1,j))*u(i,j)
                  flx(i,j)  = fldlo(i-1,j)*uloc(i,j)
                  ucumdt(i,j) = ucumdt(i,j) +  &
                                min(dt2-ucumdt(i,j),dtloc(i-1,j))
                else !u(i,j).lt.0
                  uloc(i,j) = min(dt2-ucumdt(i,j),dtloc(i,j))*u(i,j)
                  flx(i,j)  = fldlo(i,j)*uloc(i,j)
                  ucumdt(i,j) = ucumdt(i,j) +  &
                                min(dt2-ucumdt(i,j),dtloc(i,j))
                endif !u
                flxcum(i,j) = flxcum(i,j) + flx(i,j)
              else !ucumdt==dt2
                uloc(i,j) = 0.0
                flx(i,j)  = 0.0
              endif !ucumdt
            endif !iu
            if (SEA_V) then
              if (vcumdt(i,j).ne.dt2) then
                if (v(i,j).ge.0) then
                  vloc(i,j) = min(dt2-vcumdt(i,j),dtloc(i,j-1))*v(i,j)
                  fly(i,j)  = fldlo(i,j-1)*vloc(i,j)
                  vcumdt(i,j) = vcumdt(i,j) +  &
                                min(dt2-vcumdt(i,j),dtloc(i,j-1))
                else !v.lt.0
                  vloc(i,j) = min(dt2-vcumdt(i,j),dtloc(i,j))*v(i,j)
                  fly(i,j)  = fldlo(i,j)*vloc(i,j)
                  vcumdt(i,j) = vcumdt(i,j) +  &
                                min(dt2-vcumdt(i,j),dtloc(i,j))
                endif
                flycum(i,j) = flycum(i,j) + fly(i,j)
              else !vcumdt==dt2
                vloc(i,j) = 0.0
                fly(i,j)  = 0.0
             endif !vcumdt
           endif !iv
         enddo !i
        enddo !j
!$OMP   END PARALLEL DO
!
        margin = mbdy_a - 2

!$OMP   PARALLEL DO PRIVATE(j,i,qp) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-margin, jj+margin
          do i=1-margin,ii+margin
            if (SEA_P) then
              if (lcalc(i,j)) then  !lcalc from prevous iteration
                qp = hloc(i,j) - (uloc(i+1,j)-uloc(i,j)+ &
                                  vloc(i,j+1)-vloc(i,j))*scali(i,j)
!
!                                it may happen that the cfl is violated,
!                                                leading to a negative h
!                                                -----------------------
                if (qp.gt.0.0) then
                  fldlo(i,j) = ((epsil+hloc(i,j))*fldlo(i,j) -  &
                            (flx(i+1,j)-flx(i,j)+ &
                             fly(i,j+1)-fly(i,j))*scali(i,j))/(epsil+qp)
                endif !qp
                hloc(i,j)  = qp
                lcalc(i,j) = ucumdt(i+1,j).ne.dt2.or. &
                             ucumdt(i,j)  .ne.dt2.or. &
                             vcumdt(i,j+1).ne.dt2.or. &
                             vcumdt(i,j)  .ne.dt2
              endif !lcalc
            endif !ip
          enddo !i
        enddo !j
!$OMP  END PARALLEL DO
!
        call xctilr( hloc(1-nbdy,1-nbdy),1,1, mbdy_a,mbdy_a, halo_ps)
        call xctilr(fldlo(1-nbdy,1-nbdy),1,1, mbdy_a,mbdy_a, halo_ps)
!
        if     (lpipe .and. lpipe_advem) then
! ---     compare two model runs.
          call pipe_compare_sym1(fldlo,  ip,'aditer:fldlo')
          call pipe_compare_sym1(hloc,   ip,'aditer:hloc ')
          call pipe_compare_sym1(flxcum, iu,'aditer:xcum ')
          call pipe_compare_sym1(flycum, iv,'aditer:ycum ')
        endif
!
      enddo !iter
!
!.....Leapfrog step using high order scheme
!
! --- rhs: u, v, fld+, flx, fly
! --- lhs: fax, fay
!
      margin = mbdy_a - 2
!
!$OMP PARALLEL DO PRIVATE(j,i,fhx,fhy) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin, jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            fhx=u(i,j)*0.5*(fldc(i,j)+fldc(i-1,j))  ! 2nd order in space
            fax(i,j)= fhx-flxcum(i,j)/dt2           ! anti-diffusion x-flux
          endif !iu
          if (SEA_V) then
            fhy=v(i,j)*0.5*(fldc(i,j)+fldc(i,j-1))  ! 2nd order in space
            fay(i,j)= fhy-flycum(i,j)/dt2           ! anti-diffusion y-flux
          endif !iv
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!              
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.           
        call pipe_compare_sym2(fax, iu,'ad:xcum:fax ', &
                               fay, iv,'ad:ycum:fay ')
      endif                                   
!
      do j=1-margin,jj+margin
        do l=1,isp(j) !ok
          if     (ifp(j,l).ge. 1-margin) then
            fax(ifp(j,l)  ,j)=0.0
          endif
          if     (ilp(j,l).lt.ii+margin) then
            fax(ilp(j,l)+1,j)=0.0
          endif
        enddo !l
      enddo !j
!
      do i=1-margin,ii+margin
        do l=1,jsp(i) !ok
          if     (jfp(i,l).ge. 1-margin) then
            fay(i,jfp(i,l)  )=0.0
          endif
          if     (jlp(i,l).lt.jj+margin) then
            fay(i,jlp(i,l)+1)=0.0
          endif
        enddo !l
      enddo !i
!              
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.           
        call pipe_compare_sym2(fax, iu,'ad:ip:0:fax ', &
                               fay, iv,'ad:ip:0:fay ')
      endif                                   
!========================================================
!
! --- finish computation of antidiffusive fluxes
!
! --- rhs: fax+,fay+,fldlo,fcn,scal
! --- lhs: rp, rm
!
      margin = mbdy_a - 3
!
      qdt2 = 1.0/dt2
!$OMP PARALLEL DO PRIVATE(j,i,ia,ib,ja,jb, &
!$OMP                     fqmax,fqmin,famax,famin,qp,qm) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin, jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
            ia=i-1; if (ip(ia,j).eq.0) ia=i
            ib=i+1; if (ip(ib,j).eq.0) ib=i
            ja=j-1; if (ip(i,ja).eq.0) ja=j
            jb=j+1; if (ip(i,jb).eq.0) jb=j
            fqmax = max(fldlo(i,j),fldlo(ia,j ),fldlo(ib,j ), &
                                   fldlo(i, ja),fldlo(i, jb))
            fqmin = min(fldlo(i,j),fldlo(ia,j ),fldlo(ib,j ), &
                                   fldlo(i, ja),fldlo(i, jb))
            famax = max(0.0,fax(i,  j)  ) - min(0.0,fax(i+1,j)  ) + &
                    max(0.0,fay(i,  j)  ) - min(0.0,fay(i,  j+1))
            famin = max(0.0,fax(i+1,j)  ) - min(0.0,fax(i,  j)  ) + &
                    max(0.0,fay(i,  j+1)) - min(0.0,fay(i,  j)  )
            if (famax > epsil) then
              qp = (fqmax-fldlo(i,j)) *hloc(i,j)*scal(i,j)*qdt2
              if (qp < famax) then
                rp(i,j) = qp/famax  !<1.0
              else
                rp(i,j) = 1.0
              endif
            else
              rp(i,j) = 0.0
            endif
            if (famin > epsil) then
              qm = (fldlo(i,j)-fqmin) *hloc(i,j)*scal(i,j)*qdt2
              if (qm < famin) then
                rm(i,j) = qm/famin  !<1.0
              else
                rm(i,j) = 1.0
              endif
            else
              rm(i,j) = 0.0
            endif
            fmx(i,j) = fqmax  !less restrictive maximum
            fmn(i,j) = fqmin  !less restrictive minimum
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
! --- rhs: rp+, rm+
! --- lhs: fax, fay
!
      margin = mbdy_a - 4
!
!$OMP PARALLEL DO PRIVATE(j,i,fact) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin, jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            if (fax(i,j) < 0.0) then
              fact = min(rp(i-1,j),rm(i,j))
            else
              fact = min(rp(i,j),rm(i-1,j))
            endif
            fax(i,j) = fact*fax(i,j)
          endif !iu
          if (SEA_V) then
            if (fay(i,j) < 0.0) then
              fact = min(rp(i,j-1),rm(i,j))
            else
              fact = min(rp(i,j),rm(i,j-1))
            endif
            fay(i,j) = fact*fay(i,j)
          endif !iv
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(fax, iu,'ad:fact:fax ', &
                               fay, iv,'ad:fact:fay ')
      endif
!
      margin = mbdy_a - 5
!
!$OMP PARALLEL DO PRIVATE(j,i) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin, jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
!...........apply antidiffusive flux correction
            flxdiv(i,j)=((fax(i+1,j)-fax(i,j))+ &
                         (fay(i,j+1)-fay(i,j)) )*dt2*scali(i,j)
            if (hloc(i,j).gt.0.) then
               fld(i,j)=((epsil+hloc(i,j))*fldlo(i,j)-flxdiv(i,j))/ &
                         (epsil+hloc(i,j))
            else
               fld(i,j)=fldlo(i,j)
            endif !hloc
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym1(flxdiv, ip,'ad:t2c:flxdv')
        call pipe_compare_sym1(fld,    ip,'ad:fct2c:fld')
      endif
!
      return
      end subroutine advem_fct2c

      subroutine advem_fct4(fld,fldc,u,v,fco,fcn,scal,scali,dt2)
      implicit none
!
      real,    intent(in)    :: dt2
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(inout) :: fld
      real, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
               intent(in)    :: fldc,u,v,scal,scali,fco,fcn
!
! Leapfrog 4th order FCT
! S.T. Zalesak (1979): Fully multidimensional flux-corrected
! transport algorithms for fluids, J.Comput.Phys. 31 335-362.
!
! time steps are "old", "center" and "new".
!
!  fld   - scalar field, need not be >0, old input but new output
!  fldc  - scalar field at center time step
!  u,v   - mass fluxes satisfying continuity equation (old to new)
!  fco   - thickness of the layer at old    time step
!  fcn   - thickness of the layer at new    time step
!  scal  - spatial increments (squared)
!  scali - inverse of scal
!  dt2   - temporal increment (from old to new)
!
!  on return, fld's valid halo will be 0 wide.
!
      real,    parameter :: onemu=9806.e-12 !very small layer thickness
!
      real,    parameter :: ft14= 7.0/12.0,   & !4th centered inner coeff
                            ft24=-1.0/12.0   !4th centered outer coeff
!
      real    flxdn,flxdp,flydn,flydp,q
      integer i,j,l,ia,ib,ja,jb,mbdy_a,margin
!
      real :: fhx, fhy, fqmax, fqmin, famax, famin
      real :: qdt2, qp, qm, fact

      mbdy_a = mbdy_advtyp(4)  ! = 5
!
! --- compute low-order and part of antidiffusive fluxes
!
! --- rhs: u, v, fld+
! --- lhs: flx, fly, fmx, fmn
!
      margin = mbdy_a - 1
!
!$OMP PARALLEL DO PRIVATE(j,i,q,fhx,fhy,jb,ja,ib,ia) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            if (u(i,j).ge.0.0) then
              q=fld(i-1,j)
            else
              q=fld(i  ,j)
            endif
            flx(i,j)=u(i,j)*q
          endif !iu
          if (SEA_V) then
            if (v(i,j).ge.0.0) then
              q=fld(i,j-1)
            else
              q=fld(i,j  )
            endif
            fly(i,j)=v(i,j)*q
          endif !iv
          if (SEA_P) then
            ia=ipim1(i,j) !i-1 if sea and i if land
            ib=ipip1(i,j) !i+1 if sea and i if land
            ja=ipjm1(i,j) !j-1 if sea and j if land
            jb=ipjp1(i,j) !j+1 if sea and j if land
            fmx(i,j)=max(fld(i,j),fld(ia,j),fld(ib,j), &
                                  fld(i,ja),fld(i,jb))
            fmn(i,j)=min(fld(i,j),fld(ia,j),fld(ib,j), &
                                  fld(i,ja),fld(i,jb))
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
!       call pipe_compare_sym2(tx1, iu,'ad:11:tx1   ',
!    &                         ty1, iv,'ad:11:ty1   ')
!       call pipe_compare_sym2(flx, iu,'ad:11:flx   ',
!    &                         fly, iv,'ad:11:fly   ')
      endif
!
      do j=1-margin,jj+margin
        do l=1,isp(j) !ok
          if     (ifp(j,l).ge. 1-margin) then
            flx(ifp(j,l)  ,j)=0.0
          endif
          if     (ilp(j,l).lt.ii+margin) then
            flx(ilp(j,l)+1,j)=0.0
          endif
        enddo !l
      enddo !j
!
      do i=1-margin,ii+margin
        do l=1,jsp(i) !ok
          if     (jfp(i,l).ge. 1-margin) then
            fly(i,jfp(i,l)  )=0.0
          endif
          if     (jlp(i,l).lt.jj+margin) then
            fly(i,jlp(i,l)+1)=0.0
          endif
        enddo !l
      enddo !i
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(flx, iu,'ad:22:flx   ', &
                               fly, iv,'ad:33:fly   ')
      endif
!diag if     (itests.gt.0 .and. jtests.gt.0) then
!diag   i=itests
!diag   j=jtests
!diag   write (lp,'(a,2i5,f22.3/1pe39.2/0pf21.3,1pe9.2,0pf9.3, &
!diag              1pe9.2,0pf9.3/1pe39.2/0pf39.3)') &
!diag     'advem (1)',i+i0,j+j0, &
!diag     fldc(i-1,j),u(i,j),fldc(i,j-1),v(i,j), &
!diag     fldc(i,j),v(i,j+1),fldc(i,j+1),u(i+1,j),fldc(i+1,j)
!diag endif
!
! --- rhs: flx+, fly+, fld, fmx, fmn, fco, fcn, flxdiv
! --- lhs: fldlo, fmxlo, fmnlo, flxdiv
!
      margin = mbdy_a - 2
!
!$OMP PARALLEL DO PRIVATE(j,i,q) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
!...........lo order Donor Cell step
            flxdiv(i,j)=((flx(i+1,j)-flx(i,j))+ &
                         (fly(i,j+1)-fly(i,j)) )*dt2*scali(i,j)
            q=fld(i,j)*(fco(i,j)+onemu)-flxdiv(i,j)
            !max,min should only be active for very thin layers
            fldlo(i,j)=max( fmn(i,j), min( fmx(i,j), &
                                           q/(fcn(i,j)+onemu) ))
            fmxlo(i,j) = max(fld(i,j),fldc(i,j),fldlo(i,j))
            fmnlo(i,j) = min(fld(i,j),fldc(i,j),fldlo(i,j))
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym1(fldlo,  ip,'ad:610:fldlo')
        call pipe_compare_sym1(flxdiv, ip,'ad:610:flxdv')
      endif
!
!.....Leapfrog step using high order scheme
!
! --- rhs: u, v, fldi++, flx, fly
! --- lhs: fax, fay
!
      margin = mbdy_a - 2
!
!$OMP PARALLEL DO PRIVATE(j,i,q,jb,ja,ib,ia,fhx,fhy) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            if     (iu(i-1,j).eq.0 .or. iu(i+1,j).eq.0) then
              ! 2nd order time centered
              fhx=u(i,j)*0.5*(fldc(i,j)+fldc(i-1,j)) 
            else
              ! 4th order time centered
              fhx=u(i,j)*(ft14*(fldc(i,  j)+fldc(i-1,j))+ &
                          ft24*(fldc(i+1,j)+fldc(i-2,j)) )
            endif
            fax(i,j)= fhx-flx(i,j)  ! anti-diffusion x-flux
          endif !iu
          if (SEA_V) then
            if     (iv(i,j-1).eq.0 .or. iv(i,j+1).eq.0) then
              ! 2nd order time centered
              fhy=v(i,j)*0.5*(fldc(i,j)+fldc(i,j-1))
            else
              ! 4th order time centered
              fhy=v(i,j)*(ft14*(fldc(i,j)  +fldc(i,j-1))+ &
                          ft24*(fldc(i,j+1)+fldc(i,j-2)) )
            endif
            fay(i,j)= fhy-fly(i,j)  ! anti-diffusion y-flux
          endif !iv
        enddo !i
      enddo !j
!$OMP END PARALLEL DO

      do j=1-margin,jj+margin
        do l=1,isp(j) !ok
          if     (ifp(j,l).ge. 1-margin) then
            fax(ifp(j,l)  ,j)=0.0
          endif
          if     (ilp(j,l).lt.ii+margin) then
            fax(ilp(j,l)+1,j)=0.0
          endif
        enddo !l
      enddo !j
!
      do i=1-margin,ii+margin
        do l=1,jsp(i) !ok
          if     (jfp(i,l).ge. 1-margin) then
            fay(i,jfp(i,l)  )=0.0
          endif
          if     (jlp(i,l).lt.jj+margin) then
            fay(i,jlp(i,l)+1)=0.0
          endif
        enddo !l
      enddo !i
!========================================================
!
! --- finish computation of antidiffusive fluxes
!
! --- rhs: fmnlo+,fmxlo+,fax+,fay+,fldlo,fcn,scal
! --- lhs: rp, rm
!
      margin = mbdy_a - 3
!
      qdt2 = 1.0/dt2
!$OMP PARALLEL DO PRIVATE(j,i,ia,ib,ja,jb, &
!$OMP                     fqmax,fqmin,famax,famin,qp,qm) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
            ia=ipim1(i,j) !i-1 if sea and i if land
            ib=ipip1(i,j) !i+1 if sea and i if land
            ja=ipjm1(i,j) !j-1 if sea and j if land
            jb=ipjp1(i,j) !j+1 if sea and j if land
            fqmax = max(fmxlo(i,j),fmxlo(ia,j),fmxlo(ib,j), &
                                   fmxlo(i,ja),fmxlo(i,jb))
            fqmin = min(fmnlo(i,j),fmnlo(ia,j),fmnlo(ib,j), &
                                   fmnlo(i,ja),fmnlo(i,jb))
            famax = max(0.0,fax(i, j) ) - min(0.0,fax(ib,j) ) + &
                    max(0.0,fay(i, j) ) - min(0.0,fay(i, jb))
            famin = max(0.0,fax(ib,j) ) - min(0.0,fax(i, j) ) + &
                    max(0.0,fay(i, jb)) - min(0.0,fay(i, j) )
            if (famax > 0.0) then
              qp = (fqmax-fldlo(i,j)) *fcn(i,j)*scal(i,j)*qdt2
              if (qp < famax) then
                rp(i,j) = qp/famax  !<1.0
              else
                rp(i,j) = 1.0
              endif
            else
              rp(i,j) = 0.0
            endif
            if (famin > 0.0) then
              qm = (fldlo(i,j)-fqmin) *fcn(i,j)*scal(i,j)*qdt2
              if (qm < famin) then
                rm(i,j) = qm/famin  !<1.0
              else
                rm(i,j) = 1.0
              endif
            else
              rm(i,j) = 0.0
            endif
            fmx(i,j) = fqmax  !less restrictive maximum
            fmn(i,j) = fqmin  !less restrictive minimum
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
! --- rhs: rp+, rm+
! --- lhs: fax, fay
!
      margin = mbdy_a - 4
!
!$OMP PARALLEL DO PRIVATE(j,i,fact) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_U) then
            if (fax(i,j) < 0.0) then
              fact = min(rp(i-1,j),rm(i,j))
            else
              fact = min(rp(i,j),rm(i-1,j))
            endif
            fax(i,j) = fact*fax(i,j)
          endif !iu
          if (SEA_V) then
            if (fay(i,j) < 0.0) then
              fact = min(rp(i,j-1),rm(i,j))
            else
              fact = min(rp(i,j),rm(i,j-1))
            endif
            fay(i,j) = fact*fay(i,j)
          endif !iv
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym2(flx, iu,'ad:18:flx   ', &
                               fly, iv,'ad:18:fly   ')
      endif
!
!diag i=itests
!diag j=jtests
!diag write (lp,'(''advem (2)'',2i5,f22.3/1pe39.2/0pf21.3,1pe9.2,0pf9.3, &
!diag 1pe9.2,0pf9.3/1pe39.2/0pf39.3)') i,j,fldlo(i-1,j),u(i,j),fldlo(i,j-1), &
!diag v(i,j),fldlo(i,j),v(i,j+1),fldlo(i,j+1),u(i+1,j),fldlo(i+1,j)
!
! --- rhs: fax+, fay+, fld, fcn, fmx, fmn
! --- lhs: fld
!
      margin = mbdy_a - 5
!
!$OMP PARALLEL DO PRIVATE(j,i) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
            fldao(i,j) = fld(i,j)*fco(i,j)*scal(i,j)
!
!...........apply antidiffusive flux correction
            flxdiv(i,j)=((fax(i+1,j)-fax(i,j))+ &
                         (fay(i,j+1)-fay(i,j)) )*dt2*scali(i,j)
            !max,min should only be active for very thin layers
            fld(i,j)=max( fmn(i,j), min( fmx(i,j), &
                          fldlo(i,j)-flxdiv(i,j)/(fcn(i,j)+onemu) ))
!
            fldan(i,j) = fld( i,j)*fcn(i,j)*scal(i,j)
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
      if     (lpipe .and. lpipe_advem) then
! ---   compare two model runs.
        call pipe_compare_sym1(flxdiv, ip,'ad:620:flxdv')
        call pipe_compare_sym1(fld,    ip,'ad:1620:fld ')
      endif
      return
      end subroutine advem_fct4

      subroutine tsadvc(m,n)
      use mod_cb_arrays  ! HYCOM saved arrays
      implicit none
!
      integer m,n
!
! --- -----------------------------------------------------
! --- thermodynamic variable(s): advection and diffusion.
! --- -----------------------------------------------------
! --- on entry:
! --- saln(:,:,:,n) = time step t-1
! --- saln(:,:,:,m) = time step t
!
! ---  dpo(:,:,:,n) = time step t-1
! ---  dpo(:,:,:,m) = time step t
! ---  dp( :,:,:,m) = time step t   with RA time smoothing
! ---  dp( :,:,:,n) = time step t+1
!
! --- on exit:
! --- saln(:,:,:,m) = time step t   with RA time smoothing
! --- saln(:,:,:,n) = time step t+1
! --- -----------------------------------------------------
!
      logical, parameter :: lpipe_tsadvc=.false.  !extra checking (when pipe on)
!
#if defined(RELO)
      real, save, allocatable, dimension (:,:) :: &
       sold,told,q2old,q2lold
      real, save, allocatable, dimension (:,:,:) :: &
       trold
#else
      real, save, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) :: &
       sold,told,q2old,q2lold
      real, save, dimension (1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,mxtrcr) :: &
       trold
#endif
!
      logical latemp,lath3d,ldtemp,ldth3d
      integer i,isave,j,jsave,k,ktr,l,mbdy,mdf,margin
      real sminn,smaxx,flxdiv,th3d_t, &
           factor,q,dpsold,dpsmid,dpsnew, &
           dpold,dpmid,dpnew,qdpmidn,smin,ufa,ufb,vfa,vfb
      real xmin(kdm),xmax(kdm)
      real sminny(jdm),smaxxy(jdm)
!
      character*12 text,textu,textv
!
! --- for mpdata, select posdef:
! ---   1. as a power of 2
! ---   2. so that the ratio of the standard deviation to the mean of
! ---      each field is approximately the same:
! ---        0 for -saln-,   256 for -temp-, 32 for -th3d-,
! ---        0 for -tracer-,   1 for -q2-
      real       pdzero,pdtemp,pdth3d,pdq2,eps_har
      parameter (pdzero=0.0, pdtemp=256.0, pdth3d=32.0, pdq2=1.0, &
                 eps_har=1.0e-20)
!
      real harmonc,harmonz,a,aa,b,bb
# include "stmt_fns.h"
!
! --- harmonic mean - Conservative for lateral friction
! ---               - force a and b to be non-negative for safety
      harmonz(a, b )=2.0*a*b/max((a+b),2.0*eps_har)
      harmonc(aa,bb)=harmonz(max(aa,0.0),max(bb,0.0))
#if defined(RELO)
!
      if     (.not.allocated(sold)) then
        allocate( &
                 sold(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 told(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) )
        call mem_stat_add( 2*(idm+2*nbdy)*(jdm+2*nbdy) )
                 sold = r_init
                 told = r_init
        if (mxlmy) then
          allocate( &
                  q2old(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
                 q2lold(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy) )
          call mem_stat_add( 2*(idm+2*nbdy)*(jdm+2*nbdy) )
                  q2old = r_init
                 q2lold = r_init
        endif !mxlmy
        if     (ntracr.gt.0) then
          allocate( &
                 trold(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy,ntracr) )
          call mem_stat_add( (idm+2*nbdy)*(jdm+2*nbdy)*ntracr )
                 trold = r_init
        endif !ntracr
      endif !.not.allocated
#endif
!
      ldebug_tsdif = .false.
      ldebug_advem = .false.
!
      itests = itest  !local copy
      jtests = jtest  !local copy
!
      if     (btrmas) then
        onetamas(:,:,n) = oneta(:,:,n)  !use dp rather than dp' for diffusion
        onetamas(:,:,m) = oneta(:,:,n)  !use dp rather than dp' for advection, note "n"
      else
        onetamas(:,:,n) = oneta(:,:,n)  !use dp rather than dp' for diffusion
        onetamas(:,:,m) = 1.0           !use                dp' for advection
      endif
!
      uflux(:,:) = 0.0  !sets land values
      vflux(:,:) = 0.0  !sets land values
!
      mbdy = mbdy_advtyp(abs(advtyp))  ! 2-8 depending on advection scheme
!
      if     (nbdy.lt.mbdy) then
        if     (mnproc.eq.1) then
          write(lp,'(/ a,i3,a /)') &
            'error: nbdy (dimensions.h) must be at least', &
            mbdy,' for the advection scheme indicated by advtyp'
        endif
        call xcstop('tsadvc')
               stop 'tsadvc'
      endif
!
      l = mbdy
! --- dp halo is up to date
      call xctilr(saln(1-nbdy,1-nbdy,1,1),1,2*kk, l,l, halo_ps)
      call xctilr(temp(1-nbdy,1-nbdy,1,1),1,2*kk, l,l, halo_ps)
      call xctilr(th3d(1-nbdy,1-nbdy,1,1),1,2*kk, l,l, halo_ps)
      call xctilr(uflx(1-nbdy,1-nbdy,1  ),1,  kk, l,l, halo_uv)
      call xctilr(vflx(1-nbdy,1-nbdy,1  ),1,  kk, l,l, halo_vv)
      do ktr= 1,ntracr
        call xctilr(tracer(1-nbdy,1-nbdy,1,1,ktr),1,2*kk, l,l, halo_ps)
      enddo !ktr
      if (mxlmy) then
        call xctilr(q2( 1-nbdy,1-nbdy,0,1),1,2*kk+4, l,l, halo_ps)
        call xctilr(q2l(1-nbdy,1-nbdy,0,1),1,2*kk+4, l,l, halo_ps)
      endif
!
      do 81 k=1,kk
!
! --- ---------------------------------------------------
! --- advection of thermodynamic variable(s) (and tracer)
! --- ---------------------------------------------------
!
! --- for isopycnic vertical coordinates:
! ---   advect -th3d- and -saln- in the mixed layer (layer 1),
! ---   advect            -saln- only in all other layers
! --- for hybrid vertical coordinates:
! ---   advect -temp- and -saln- in all layers if advflg==0,
! ---   advect -th3d- and -saln- in all layers if advflg==1
!
      latemp =  k.le.nhybrd .and. advflg.eq.0       ! advect temp
      lath3d = (k.le.nhybrd .and. advflg.eq.1) .or. &
               (k.eq.1      .and. isopyc     )      ! advect th3d
!
! --- smooth mixed-layer mass fluxes in lateral direction
      if(isopyc .and. k.eq.1) then
!
! ---   rhs: vflx+, uflx+
! ---   lhs: vflux, uflux
!
        margin = mbdy - 1
!
        do j=1-margin,jj+margin
          do i=1-margin,ii+margin
            if (SEA_V) then
              if (iv(i-1,j).ne.0) then
                vfa=vflx(i-1,j,1)
              else 
                vfa=vflx(i,  j,1)
              endif
              if (iv(i+1,j).ne.0) then
                vfb=vflx(i+1,j,1)
              else
                vfb=vflx(i,  j,1)
              endif
              vflux(i,j)=.5*vflx(i,j,1)+.25*(vfa+vfb)
            endif !iv
            if (SEA_U) then
              if (iu(i,j-1).ne.0) then
                ufa=uflx(i,j-1,1)
              else
                ufa=uflx(i,j,  1)
              endif
              if (iu(i,j+1).ne.0) then
                ufb=uflx(i,j+1,1)
              else
                ufb=uflx(i,j,  1)
              endif
              uflux(i,j)=.5*uflx(i,j,1)+.25*(ufa+ufb)
            endif !iu
          enddo !i
        enddo !j
      endif
!
! ---   rhs: temp, saln, uflux+, vflux+, dp
! ---   lhs: told, sold, util1, util2, util3, temp, th3d
!
! ---   util1 = fco = thickness of the layer at old    time step
! ---   util2 = fcn = thickness of the layer at new    time step
!
        margin = mbdy - 1  ! util[12] at mbdy-2
!
!$OMP PARALLEL DO PRIVATE(j,i,ktr,flxdiv) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
!
! ---       save for time smoothing
            if     (latemp) then
              told(i,j)=temp(i,j,k,n)
            elseif (lath3d) then
              told(i,j)=th3d(i,j,k,n)
            endif
            sold(i,j)=saln(i,j,k,n)
            do ktr= 1,ntracr
              trold(i,j,ktr)=tracer(i,j,k,n,ktr)
            enddo
            if (mxlmy) then
              q2old( i,j)=q2( i,j,k,n)
              q2lold(i,j)=q2l(i,j,k,n)
            endif
!
! --- before calling 'advem', make sure (a) mass fluxes are consistent
! --- with layer thickness change, and (b) all fields are positive-definite
            if(isopyc .and. k.eq.1) then
              flxdiv=((uflux(i+1,j)  -uflux(i,j)  ) &
                     +(vflux(i,j+1)  -vflux(i,j)  ))*delt1*scp2i(i,j)
            else
              flxdiv=((uflx( i+1,j,k)-uflx( i,j,k)) &
                     +(vflx( i,j+1,k)-vflx( i,j,k)))*delt1*scp2i(i,j)
            endif
            util1(i,j)=max(onetamas(i,j,m)*dp(i,j,k,n)+flxdiv,0.0)  !old, note "m"
            util2(i,j)=max(onetamas(i,j,m)*dp(i,j,k,n),       0.0)  !new, note "m"
          endif !ip
        enddo !i
      enddo !j
!$OMP END PARALLEL DO
!
      if     (lpipe .and. lpipe_tsadvc) then
! ---   compare two model runs.
        write(text,'(a10,i2)') '49:sold,k=',k
        call pipe_compare_sym1(sold,ip,text)
        write(text,'(a10,i2)') '49:told,k=',k
        call pipe_compare_sym1(told,ip,text)
        write(text,'(a10,i2)') '49:utl1,k=',k
        call pipe_compare_sym1(util1,ip,text)
        write(text,'(a10,i2)') '49:utl2,k=',k
        call pipe_compare_sym1(util2,ip,text)
        write (textu,'(a9,i3)') 'uflx   k=',k
        write (textv,'(a9,i3)') 'vflx   k=',k
        call pipe_compare_sym2(uflx(1-nbdy,1-nbdy,k),  iu,textu, &
                               vflx(1-nbdy,1-nbdy,k),  iv,textv)
        write (text,'(a9,i3)') 'temp.n k=',k
        call pipe_compare_sym1(temp(1-nbdy,1-nbdy,k,n),ip,text)
        write (text,'(a9,i3)') 'saln.n k=',k
        call pipe_compare_sym1(saln(1-nbdy,1-nbdy,k,n),ip,text)
        write (text,'(a9,i3)') 'th3d.n k=',k
        call pipe_compare_sym1(th3d(1-nbdy,1-nbdy,k,n),ip,text)
      endif
!
! --- rhs: temp.[mn], th3d.[mn], saln.[mn], uflx, vflx, util[12]
! --- lhs: temp.n, th3d.n, saln.n
!
      if     (latemp) then
        call advem(advtyp,temp( 1-nbdy,1-nbdy,k,n), &
                          temp( 1-nbdy,1-nbdy,k,m), &
                          uflx( 1-nbdy,1-nbdy,k), &
                          vflx( 1-nbdy,1-nbdy,k), &
                          util1,util2, &
                          pdtemp, scp2,scp2i,delt1, btrmas)
        call advem(advtyp,saln( 1-nbdy,1-nbdy,k,n), &
                          saln( 1-nbdy,1-nbdy,k,m), &
                          uflx( 1-nbdy,1-nbdy,k), &
                          vflx( 1-nbdy,1-nbdy,k), &
                          util1,util2, &
                          pdzero, scp2,scp2i,delt1, btrmas)
      elseif (lath3d .and. hybrid) then
        call advem(advtyp,th3d( 1-nbdy,1-nbdy,k,n), &
                          th3d( 1-nbdy,1-nbdy,k,m), &
                          uflx( 1-nbdy,1-nbdy,k), &
                          vflx( 1-nbdy,1-nbdy,k), &
                          util1,util2, &
                          pdth3d, scp2,scp2i,delt1, btrmas)
        call advem(advtyp,saln( 1-nbdy,1-nbdy,k,n), &
                          saln( 1-nbdy,1-nbdy,k,m), &
                          uflx( 1-nbdy,1-nbdy,k), &
                          vflx( 1-nbdy,1-nbdy,k), &
                          util1,util2, &
                          pdzero, scp2,scp2i,delt1, btrmas)
      elseif (lath3d .and. isopyc) then  ! MICOM-like upper layer
        call advem(advtyp,th3d( 1-nbdy,1-nbdy,k,n), &
                          th3d( 1-nbdy,1-nbdy,k,m), &
                          uflux, &
                          vflux, &
                          util1,util2, &
                          pdth3d, scp2,scp2i,delt1, btrmas)
        call advem(advtyp,saln( 1-nbdy,1-nbdy,k,n), &
                          saln( 1-nbdy,1-nbdy,k,m), &
                          uflux, &
                          vflux, &
                          util1,util2, &
                          pdzero, scp2,scp2i,delt1, btrmas)
      else   ! exactly isopycnal layer
        call advem(advtyp,saln( 1-nbdy,1-nbdy,k,n), &
                          saln( 1-nbdy,1-nbdy,k,m), &
                          uflx( 1-nbdy,1-nbdy,k), &
                          vflx( 1-nbdy,1-nbdy,k), &
                          util1,util2, &
                          pdzero, scp2,scp2i,delt1, btrmas)
      endif
      do ktr= 1,ntracr
        if     (trcflg(ktr).eq.2) then !temperature tracer
          call advem(advtyp,tracer(1-nbdy,1-nbdy,k,n,ktr), &
                            tracer(1-nbdy,1-nbdy,k,m,ktr), &
                            uflx(  1-nbdy,1-nbdy,k), &
                            vflx(  1-nbdy,1-nbdy,k), &
                            util1,util2, &
                            pdtemp, scp2,scp2i,delt1, btrmas)
        else
!         ldebug_advem = k.eq.kk !debug layer kk tracer advection
          call advem(advtyp,tracer(1-nbdy,1-nbdy,k,n,ktr), &
                            tracer(1-nbdy,1-nbdy,k,m,ktr), &
                            uflx(  1-nbdy,1-nbdy,k), &
                            vflx(  1-nbdy,1-nbdy,k), &
                            util1,util2, &
                            pdzero, scp2,scp2i,delt1, btrmas)
!         ldebug_advem = .false. !turn off debugging of advection
        endif !trcflg
      enddo !ktr
      if (mxlmy) then
        call advem(advtyp,q2(   1-nbdy,1-nbdy,k,n), &
                          q2(   1-nbdy,1-nbdy,k,m), &
                          uflx( 1-nbdy,1-nbdy,k), &
                          vflx( 1-nbdy,1-nbdy,k), &
                          util1,util2, &
                          pdq2,   scp2,scp2i,delt1, btrmas)
        call advem(advtyp,q2l(  1-nbdy,1-nbdy,k,n), &
                          q2l(  1-nbdy,1-nbdy,k,m), &
                          uflx( 1-nbdy,1-nbdy,k), &
                          vflx( 1-nbdy,1-nbdy,k), &
                          util1,util2, &
                          pdq2,   scp2,scp2i,delt1, btrmas)
      endif
!
      if     (lpipe .and. lpipe_tsadvc) then
! ---   compare two model runs.
        write (text,'(a9,i3)') 'temp.n k=',k
        call pipe_compare_sym1(temp(1-nbdy,1-nbdy,k,n),ip,text)
        write (text,'(a9,i3)') 'saln.n k=',k
        call pipe_compare_sym1(saln(1-nbdy,1-nbdy,k,n),ip,text)
        write (text,'(a9,i3)') 'th3d.n k=',k
        call pipe_compare_sym1(th3d(1-nbdy,1-nbdy,k,n),ip,text)
      endif
!
!diag if (itests.gt.0.0and.jtests.gt.0) &
!diag write (lp,'(i9,2i5,i3,'' th,s,dp after advection  '',2f9.3,f8.2)') &
!diag nstep,itests,jtests,k,temp(itests,jtests,k,n),saln(itests,jtests,k,n), &
!diag dp(itests,jtests,k,n)*qonem
!
      if     (mod(nstep,3).eq.0 .or. diagno) then
!$OMP   PARALLEL DO PRIVATE(j,i) &
!$OMP            SCHEDULE(STATIC,jblk) !NOCSD
        do j=1,jj
          sminny(j)= 999.  !simplifies OpenMP parallelization
          smaxxy(j)=-999.  !simplifies OpenMP parallelization
!DIR$     PREFERVECTOR
          do i=1,ii
            if (SEA_P) then
              if     (dp(i,j,k,n).gt.onemm) then
                sminny(j)=min(sminny(j),saln(i,j,k,n))
                smaxxy(j)=max(smaxxy(j),saln(i,j,k,n))
              endif
            endif !ip
          enddo !i
        enddo !j
!$OMP   END PARALLEL DO !NOCSD
        xmin(k) = minval(sminny(1:jj))
        xmax(k) = maxval(smaxxy(1:jj))
      endif !every 3 time steps or diagno
!
 81   continue  ! k=1,kk
!
      call pipe_comparall(m,n, 'advem,  step')
!
! --- check for negative scalar fields.
!
      if     (mod(nstep,3).eq.0 .or. diagno) then
        call xcminr(xmin(1:kk))
        call xcmaxr(xmax(1:kk))
!
        do k= 1,kk
          sminn=xmin(k)
          smaxx=xmax(k)
!
          if (sminn.lt.0.0) then
            do j=1,jj
              do i=1,ii
                if (SEA_P) then
                  if (saln(i,j,k,n).eq.sminn) then
                    write (lp,'(i9,a,2i5,i3,a,f10.2)')  &
                      nstep,' i,j,k =',i+i0,j+j0,k, &
                      ' neg. saln after advem call ', &
                      saln(i,j,k,n)
                  endif !sminn
                endif !ip
              enddo !i
            enddo !j
            call xcsync(flush_lp)
          endif
!
          if (diagno) then
            if     (mnproc.eq.1) then
            if     (sminn.le.smaxx) then
              write (lp,'(i9,i3, a,2f9.3, a,1pe9.2,a)') &
                nstep,k, &
                ' min/max of s after advection:',sminn,smaxx, &
                '   (range:',smaxx-sminn,')'
            else
              write (lp,'(i9,i3, a,a)') &
                nstep,k, &
                ' min/max of s after advection:',' N/A (thin layer)'
            endif !normal:thin layer
            call flush(lp)
            endif
          endif
        enddo !k
      endif !every 3 time steps or diagno
!
! --- --------------------------------------
! --- diffusion of thermodynamic variable(s)
! --- --------------------------------------
!
      if     (temdf2.gt.0.0) then
        mdf = 2  !Laplacian
        call xctilr(saln(   1-nbdy,1-nbdy,1,n),1,kk, mdf,mdf, halo_ps)
        call xctilr(temp(   1-nbdy,1-nbdy,1,n),1,kk, mdf,mdf, halo_ps)
        call xctilr(th3d(   1-nbdy,1-nbdy,1,n),1,kk, mdf,mdf, halo_ps)
        if (mxlmy) then
          call xctilr(q2( 1-nbdy,1-nbdy,0,n),1,kk+2, mdf,mdf, halo_ps)
          call xctilr(q2l(1-nbdy,1-nbdy,0,n),1,kk+2, mdf,mdf, halo_ps)
        endif
        do ktr= 1,ntracr
          call xctilr(tracer(1-nbdy,1-nbdy,1,n,ktr), &
                                               1,kk, mdf,mdf, halo_ps)
        enddo !ktr
!
        do k=1,kk
!
! ---     for isopycnic vertical coordinates:
! ---       diffuse -th3d- and -saln- in the mixed layer (layer 1),
! ---       diffuse            -saln- only in all other layers
! ---     for hybrid vertical coordinates:
! ---       diffuse -saln- in all layers
! ---       diffuse -temp- in all layers if temdfc < 0.0
! ---       diffuse -th3d- in all layers if temdfc < 1.0
! ---       if 0.0 < temdfc < 1.0:
! ---         combine -temp- and -th3d- diffusion in density space
!
          ldtemp =  k.le.nhybrd .and. temdfc.gt.0.0     ! diffus temp
          ldth3d = (k.le.nhybrd .and. temdfc.lt.1.0) .or. &
                   (k.eq.1      .and. isopyc     )      ! diffus th3d
          if     (ldtemp .and. ldth3d) then ! diffus temp and th3d
            call tsdff_2x(th3d(1-nbdy,1-nbdy,k,n), &
                          temp(1-nbdy,1-nbdy,k,n))
            call tsdff_1x(saln(1-nbdy,1-nbdy,k,n))
          elseif (ldtemp) then ! diffus temp
            call tsdff_2x(temp(1-nbdy,1-nbdy,k,n), &
                          saln(1-nbdy,1-nbdy,k,n))
          elseif (ldth3d) then ! diffus th3d
            call tsdff_2x(th3d(1-nbdy,1-nbdy,k,n), &
                          saln(1-nbdy,1-nbdy,k,n))
          else   ! exactly isopycnal layer
            call tsdff_1x(saln(1-nbdy,1-nbdy,k,n))
          endif
          if (mxlmy) then
            call tsdff_2x(q2(  1-nbdy,1-nbdy,k,n), &
                          q2l( 1-nbdy,1-nbdy,k,n))
          endif !mxlmy
          do ktr= 1,ntracr,2
            if     (ktr+1.le.ntracr) then
              call tsdff_2x(tracer(1-nbdy,1-nbdy,k,n,ktr), &
                            tracer(1-nbdy,1-nbdy,k,n,ktr+1))
            else
!             ldebug_tsdif = k.eq.kk !debug layer kk tracer diffusion
              call tsdff_1x(tracer(1-nbdy,1-nbdy,k,n,ktr))
!             ldebug_tsdif = .false. !turn off debugging of diffusion
            endif
          enddo !ktr
        enddo !k
!
!       non-independent thermodynamic variable
!
        margin = 0
!$OMP   PARALLEL DO PRIVATE(j,k,i,ldtemp,ldth3d,th3d_t) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-margin,jj+margin
          do k=1,kk
            ldtemp =  k.le.nhybrd .and. temdfc.gt.0.0     ! diffus temp
            ldth3d = (k.le.nhybrd .and. temdfc.lt.1.0) .or. &
                     (k.eq.1      .and. isopyc     )      ! diffus th3d
            do i=1-margin,ii+margin
              if (SEA_P) then
                if     (ldtemp .and. ldth3d) then
! ---             combine -temp- and -th3d- diffusion in density space
                  th3d_t       =sig(temp(i,j,k,n),saln(i,j,k,n))-thbase
                  th3d(i,j,k,n)=(1.0-temdfc)*th3d(i,j,k,n) + &
                                     temdfc *th3d_t
                  temp(i,j,k,n)=tofsig(th3d(i,j,k,n)+thbase, &
                                       saln(i,j,k,n))
                elseif (ldtemp) then
                  th3d(i,j,k,n)=sig(temp(i,j,k,n),saln(i,j,k,n))-thbase
                elseif (ldth3d) then
                  temp(i,j,k,n)=tofsig(th3d(i,j,k,n)+thbase, &
                                       saln(i,j,k,n))
                else   ! exactly isopycnal layer
                  th3d(i,j,k,n)=theta(i,j,k)
                  temp(i,j,k,n)=tofsig(th3d(i,j,k,n)+thbase, &
                                       saln(i,j,k,n))
                endif !ld....
              endif !ip
            enddo !i
          enddo !k
        enddo !j
!$OMP   END PARALLEL DO
      endif !temdf2.gt.0.0
!
      do k=1,kk
        if     (lpipe .and. lpipe_tsadvc) then
! ---     compare two model runs.
          write (text,'(a9,i3)') 'util1  k=',k
          call pipe_compare_sym1(util1,ip,text)
          write (text,'(a9,i3)') 'util2  k=',k
          call pipe_compare_sym1(util2,ip,text)
          write (text,'(a9,i3)') 'temp.n k=',k
          call pipe_compare_sym1(temp(1-nbdy,1-nbdy,k,n),ip,text)
          write (text,'(a9,i3)') 'saln.n k=',k
          call pipe_compare_sym1(saln(1-nbdy,1-nbdy,k,n),ip,text)
          write (text,'(a9,i3)') 'th3d.n k=',k
          call pipe_compare_sym1(th3d(1-nbdy,1-nbdy,k,n),ip,text)
        endif
!
!diag   if (itests.gt.0.and.jtests.gt.0) then &
!diag     write (lp,'(i9,2i5,i3,a,2f9.3,f8.2)') &
!diag       nstep,itests+i0,jtests+j0,k, &
!diag       ' t,s,dp after isopyc.mix.', &
!diag       temp(itests,jtests,k,n),saln(itests,jtests,k,n), &
!diag       dp(itests,jtests,k,n)*qonem &
!diag     call flush(lp) &
!diag   endif
!
      enddo !k

      return
 
      contains
!
      subroutine tsdff_2x(fld1,fld2)
      real fld1(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy), &
           fld2(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)
!
! --- Laplacian diffusion for two scalar fields
!
        margin = 1
!
!$OMP   PARALLEL DO PRIVATE(j,i,factor) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-margin,jj+margin
          do i=1-margin,ii+margin
            if (SEA_U) then
              factor=temdf2*aspux(i,j)* &
                     scuy(i,j)*harmonc(dp(i-1,j,k,n)*onetamas(i-1,j,n) &
                                      ,dp(i  ,j,k,n)*onetamas(i,  j,n))
              uflux (i,j)=factor*(fld1(i-1,j)-fld1(i,j))
              uflux2(i,j)=factor*(fld2(i-1,j)-fld2(i,j))
            endif !iu
            if (SEA_V) then
              factor=temdf2*aspvy(i,j)* &
                     scvx(i,j)*harmonc(dp(i,j-1,k,n)*onetamas(i,j-1,n) &
                                      ,dp(i,j  ,k,n)*onetamas(i,j  ,n))
              vflux (i,j)=factor*(fld1(i,j-1)-fld1(i,j))
              vflux2(i,j)=factor*(fld2(i,j-1)-fld2(i,j))
            endif !iv
          enddo !i
        enddo !j
!$OMP   END PARALLEL DO
!
        if     (lpipe .and. lpipe_tsadvc) then
! ---     compare two model runs.
          write (textu,'(a9,i3)') 'uflux  k=',k
          write (textv,'(a9,i3)') 'vflux  k=',k
          call pipe_compare_sym2(uflux, iu,textu, &
                                 vflux, iv,textv)
          write (textu,'(a9,i3)') 'uflux2 k=',k
          write (textv,'(a9,i3)') 'vflux2 k=',k
          call pipe_compare_sym2(uflux2,iu,textu, &
                                 vflux2,iv,textv)
        endif
!
! ---   rhs: dp.n, uflux+, vflux+, uflux2+, vflux2+
! ---   lhs: saln.n, temp.n, th3d.n
!
        margin = 0
!
!$OMP   PARALLEL DO PRIVATE(j,i,factor) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-margin,jj+margin
          do i=1-margin,ii+margin
            if (SEA_P) then
              factor=-delt1/(scp2(i,j)*max(dp(i,j,k,n)*onetamas(i,j,n), &
                                           eps_har))
              util1(i,j)=((uflux (i+1,j)-uflux (i,j)) &
                         +(vflux (i,j+1)-vflux (i,j)))*factor
              fld1(i,j)=fld1(i,j)+util1(i,j)
              util2(i,j)=((uflux2(i+1,j)-uflux2(i,j)) &
                         +(vflux2(i,j+1)-vflux2(i,j)))*factor
              fld2(i,j)=fld2(i,j)+util2(i,j)
!
!diag         if (i.eq.itests.and.j.eq.jtests) then
!diag           if (1.le.i .and. i.le.ii .and. &
!diag               1.le.j .and. j.le.jj      ) then &
!diag             write (lp,100) nstep,i+i0,j+j0,k,'t,s,dt,ds', &
!diag             fld1(i,j),fld2(i,j),util1(i,j),util2(i,j) &
!diag             call flush(lp) &
!diag           endif &
!diag         endif
!
            endif !ip
          enddo !i
        enddo !j
!$OMP   END PARALLEL DO
!
        if     (lpipe .and. lpipe_tsadvc) then
! ---     compare two model runs.
          write (text,'(a9,i3)') 'util1  k=',k
          call pipe_compare_sym1(util1,ip,text)
          write (text,'(a9,i3)') 'util2  k=',k
          call pipe_compare_sym1(util2,ip,text)
          write (text,'(a9,i3)') 'fld1.n k=',k
          call pipe_compare_sym1(fld1(1-nbdy,1-nbdy),ip,text)
          write (text,'(a9,i3)') 'fld2.n k=',k
          call pipe_compare_sym1(fld2(1-nbdy,1-nbdy),ip,text)
        endif
!
      return
      end subroutine tsdff_2x
!
      subroutine tsdff_1x(fld1)
      real fld1(1-nbdy:idm+nbdy,1-nbdy:jdm+nbdy)
!
! --- Laplacian  diffusion for a single scalar field
!
        margin = 1
!
!$OMP   PARALLEL DO PRIVATE(j,i,factor) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-margin,jj+margin
          do i=1-margin,ii+margin
            if (SEA_U) then
              factor=temdf2*aspux(i,j)* &
                     scuy(i,j)*harmonc(dp(i-1,j,k,n)*onetamas(i-1,j,n) &
                                      ,dp(i  ,j,k,n)*onetamas(i,  j,n))
              uflux (i,j)=factor*(fld1(i-1,j)-fld1(i,j))
            endif !iu
            if (SEA_V) then
              factor=temdf2*aspvy(i,j)* &
                     scvx(i,j)*harmonc(dp(i,j-1,k,n)*onetamas(i,j-1,n) &
                                      ,dp(i,j  ,k,n)*onetamas(i,j  ,n))
              vflux (i,j)=factor*(fld1(i,j-1)-fld1(i,j))
            endif !iv
          enddo !i
        enddo !j
!$OMP   END PARALLEL DO
!
        if     (lpipe .and. lpipe_tsadvc) then
! ---     compare two model runs.
          write (textu,'(a9,i3)') 'uflux  k=',k
          write (textv,'(a9,i3)') 'vflux  k=',k
          call pipe_compare_sym2(uflux, iu,textu, &
                                 vflux, iv,textv)
        endif
!
      if     (ldebug_tsdif .and. itests.gt.0 .and. jtests.gt.0) then
        i=itests
        j=jtests
        if (SEA_U) then
          write (lp,'(a,2i5,3e16.6)') &
            'tsdff: u dp ',i+i0,j+j0, &
            uflux(i,j),  dp(i-1,j,k,n)*qonem,dp(i  ,j,k,n)*qonem
          write (lp,'(a,2i5,3e16.6)') &
            'tsdff: u fld',i+i0,j+j0, &
            uflux(i,j),fld1(i-1,j    ),    fld1(i  ,j    )
        else
          write (lp,'(a,2i5,1e16.6)') &
            'tsdff: u LND',i+i0,j+j0, &
            uflux(i,j)
        endif
        if (SEA_V) then
          write (lp,'(a,2i5,3e16.6)') &
            'tsdff: v dp ',i+i0,j+j0, &
            vflux(i,j),  dp(i,j-1,k,n)*qonem,dp(i,j  ,k,n)*qonem
          write (lp,'(a,2i5,3e16.6)') &
            'tsdff: v fld',i+i0,j+j0, &
            vflux(i,j),fld1(i,j-1    ),    fld1(i,j      )
        else
          write (lp,'(a,2i5,1e16.6)') &
            'tsdff: v LND',i+i0,j+j0, &
            vflux(i,j)
        endif
        i=itests+1
        j=jtests
        if (SEA_U) then
          write (lp,'(a,2i5,3e16.6)') &
            'tsdff: u dp ',i+i0,j+j0, &
            uflux(i,j),  dp(i-1,j,k,n)*qonem,dp(i  ,j,k,n)*qonem
          write (lp,'(a,2i5,3e16.6)') &
            'tsdff: u fld',i+i0,j+j0, &
            uflux(i,j),fld1(i-1,j    ),    fld1(i  ,j    )
        else
          write (lp,'(a,2i5,1e16.6)') &
            'tsdff: u LND',i+i0,j+j0, &
            uflux(i,j)
        endif
        i=itests
        j=jtests+1
        if (SEA_V) then
          write (lp,'(a,2i5,3e16.6)') &
            'tsdff: v dp ',i+i0,j+j0, &
            vflux(i,j),  dp(i,j-1,k,n)*qonem,dp(i,j  ,k,n)*qonem
          write (lp,'(a,2i5,3e16.6)') &
            'tsdff: v fld',i+i0,j+j0, &
            vflux(i,j),fld1(i,j-1    ),    fld1(i,j      )
        else
          write (lp,'(a,2i5,1e16.6)') &
            'tsdff: u LND',i+i0,j+j0, &
            vflux(i,j)
        endif
      endif
!
! ---   rhs: dp.n, uflux+, vflux+, uflux2+, vflux2+
! ---   lhs: saln.n, temp.n, th3d.n
!
        margin = 0
!
!$OMP   PARALLEL DO PRIVATE(j,i,factor) &
!$OMP            SCHEDULE(STATIC,jblk)
        do j=1-margin,jj+margin
          do i=1-margin,ii+margin
            if (SEA_P) then
              factor=-delt1/(scp2(i,j)*max(dp(i,j,k,n)*onetamas(i,j,n), &
                                           eps_har))
              util1(i,j)=((uflux (i+1,j)-uflux (i,j)) &
                         +(vflux (i,j+1)-vflux (i,j)))*factor
              fld1(i,j)=fld1(i,j)+util1(i,j)
!
!diag         if (i.eq.itests.and.j.eq.jtests) then
!diag           if (1.le.i .and. i.le.ii .and. &
!diag               1.le.j .and. j.le.jj      ) then &
!diag             write (lp,100) nstep,i+i0,j+j0,k,'t,s,dt,ds', &
!diag             fld1(i,j),0.0,util1(i,j),0.0 &
!diag             call flush(lp) &
!diag           endif &
!diag         endif
!
            endif !ip
          enddo !i
        enddo !j
!$OMP   END PARALLEL DO
!
        if     (lpipe .and. lpipe_tsadvc) then
! ---     compare two model runs.
          write (text,'(a9,i3)') 'util1  k=',k
          call pipe_compare_sym1(util1,ip,text)
          write (text,'(a9,i3)') 'fld1.n k=',k
          call pipe_compare_sym1(fld1(1-nbdy,1-nbdy),ip,text)
        endif
!
      if     (ldebug_tsdif .and. itests.gt.0 .and. jtests.gt.0) then
        i=itests
        j=jtests
        write (lp,'(a,2i5,3e16.6)') &
          'tsdff: p fld',i+i0,j+j0, &
          fld1(i,j),util1(i,j), &
          -delt1/(scp2(i,j)*max(dp(i,j,k,n)*onetamas(i,j,n),eps_har))
      endif
!
      return
      end subroutine tsdff_1x

!-----end contains

      end subroutine tsadvc
!
      end module mod_tsadvc
!
!  Revision history:
!
!> June 1995 - eliminated setting of salinity in massless layers (loop 46)
!>             (this is now done in mxlayr.f)
!> Aug. 1995 - omitted t/s/dp time smoothin, case of abrupt mxlayr.thk.change
!> Sep. 1995 - increased temdf2 if mixed layer occupies >90% of column
!> Mar. 2000 - removed 'cushn' and added logic to assure global conservation
!> Apr. 2000 - conversion to SI units
!> Apr. 2000 - changed i/j loop nesting to j/i
!> Aug. 2000 - temp advection and diffusion only for hybrid vertical coordinate
!> Nov. 2000 - nhybrd T&S advection layers, kdm-nhybrd S advection layers
!> Nov. 2000 - T&S or th&S advection/diffusion based on advflg
!> Feb. 2001 - placed advem in a module
!> May  2002 - diffusion coefficent based on max(sc?x,sc?y)
!> Aug. 2003 - separate PCM and MPDATA versions (advtyp)
!> Aug. 2003 - added FCT2 and UTOPIA advection options.
!> Nov. 2003 - per layer diffusion routine for 1 or 2 scalar fields
!> Feb. 2008 - fixed famin,famax bugs in FCT2/4
!> Jun. 2008 - fixed q2,q2l halo update bug
!> Aug  2011 - reworked Robert-Asselin filter, RA of dp now in cnuity
!> Aug  2012 - RA filter now exactly conserves constant tracers
!> Sep  2012 - RA filter not applied if thickness < onezm
!> Apr. 2014 - replace ip with ipa for mass sums
!> May  2014 - use land/sea masks (e.g. ip) to skip land
!> May  2014 - use ipim1,ipip1,ipjm1,ipjp1 as sea-only neighbors
!> Aug. 2018 - Robert-Asselin filter now in mod_asselin
!> Aug. 2018 - btrmas added, use onetamas to simplify logic
!> Aug. 2018 - replaced itest,jtest with itests,jtests
!> Nov. 2018 - always use oneta for diffusion
!> Nov. 2018 - replaced min(1.0,A/B) with if(A<B) to avoid overflows
