#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_asselin
      use mod_xc         ! HYCOM communication interface
      use mod_cb_arrays  ! HYCOM saved arrays
      use mod_pipe       ! HYCOM debugging interface
!
      implicit none
!
! --- module for HYCOM Robert-Asselin filter of scalar fields
!
      private !! default is private
      public  :: asselin_save,asselin_filter

      contains

      subroutine asselin_save(m,n)
!
      integer m,n
!
! --- save time level t-1 for Robert-Asselin filter
! --- dpo is initialized in cnuity
!
! --- on exit:
! ---    otemp(:,:,:)   = time step t-1
! ---    osaln(:,:,:)   = time step t-1
! ---    oth3d(:,:,:)   = time step t-1
! ---  otracer(:,:,:,:) = time step t-1
!
! ---   onetao(:,:,  n) = time step t-1
! ---   oneta( :,:,  n) = time step t-1
! ---   onetao(:,:,  m) = time step t
! ---   oneta( :,:,  m) = time step t
!
      integer i,j,k,ktr
!
!$OMP PARALLEL DO PRIVATE(j,i,k,ktr)
      do j=1,jj
        do i=1,ii
          if (SEA_P) then
          oneta( i,j,n) = max( oneta0, 1.0 + pbavg(i,j,n)/pbot(i,j) )  !t-1
          oneta( i,j,m) = max( oneta0, 1.0 + pbavg(i,j,m)/pbot(i,j) )  !t
          onetao(i,j,n) = oneta( i,j,n)
          onetao(i,j,m) = oneta( i,j,m)
          endif
        enddo !i
        do k= 1,kk
          do i=1,ii
            otemp(i,j,k)   = temp(i,j,k,n)
            osaln(i,j,k)   = saln(i,j,k,n)
            oth3d(i,j,k)   = th3d(i,j,k,n)
            do ktr= 1,ntracr
              otracer(i,j,k,ktr) = tracer(i,j,k,n,ktr)
            enddo !ktr
          enddo !i
        enddo !k
        if     (mxlmy) then
          do k= 1,kk
            do i=1,ii
               oq2(i,j,k) =  q2(i,j,k,n)
              oq2l(i,j,k) = q2l(i,j,k,n)
            enddo !i
          enddo !k
        endif !mxkmy
      enddo !j
!
      call xctilr(oneta( 1-nbdy,1-nbdy,1),1,2, 6,6, halo_ps)
      call xctilr(onetao(1-nbdy,1-nbdy,1),1,2, 6,6, halo_ps)
!
      return
      end subroutine asselin_save

      subroutine asselin_filter(m,n)
!
      integer m,n
!
! --- Apply Robert-Asselin filter
!
      logical, parameter ::  lpipe_asselin=.false.  !extra checking (when pipe on)
      logical, parameter :: ldebug_asselin=.false.  !debugging, usually false
!
      real,    parameter :: onezm=9806.e-20 !zepto-m, insignificant thickness
!
      logical          latemp,lath3d
      integer          i,j,k,ktr,margin
      real             q,dpsold,dpsmid,dpsnew, &
                       dpold,dpmid,dpmidn,dpnew,qdpmidn,smin
      double precision ssum(4),s1(4)
!
      character*12 text
!
# include "stmt_fns.h"
!
      margin = 0  !at end of time step
!
      if     (ldebug_asselin) then
        ssum(1:4) = 0.0d0
      endif !debug
!
!$OMP PARALLEL DO PRIVATE(j,k,i,ktr,dpsold,dpsmid,dpsnew,q, &
!$OMP                     dpold,dpmid,dpmidn,dpnew,qdpmidn,smin) &
!$OMP          SCHEDULE(STATIC,jblk)
      do j=1-margin,jj+margin
        do i=1-margin,ii+margin
          if (SEA_P) then
            oneta(i,j,n) = max( oneta0, 1.0 + pbavg(i,j,n)/pbot(i,j) )  !t+1
            oneta(i,j,m) = max( oneta0, 1.0 + pbavg(i,j,m)/pbot(i,j) )  !t with RA
          endif !ip
        enddo !i
        do k= 1,kk
          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
          do i=1-margin,ii+margin
            if (SEA_P) then
!
! ---         Robert-Asselin time filter of scalar fields
! ---         Note that this is smoothing oneta * dp * scalar
!
              dpold  = dpo(i,j,k,n)*onetao(i,j,n)  !t-1
              dpmid  = dpo(i,j,k,m)*onetao(i,j,m)  !t
              dpnew  = dp( i,j,k,n)*oneta( i,j,n)  !t+1
              q      = 0.5*ra2fac*(dpold+dpnew-2.0*dpmid)
              dpmidn = dpmid + q
              dp(i,j,k,m) = dpmidn/oneta(i,j,m)  !t with RA
              if     (dpmidn.gt.onezm) then !effectively non-zero
                if     (ldebug_asselin .and. i.eq.itest &
                                       .and. j.eq.jtest) then
                  ssum(1) = ssum(1) + dpold*osaln(i,j,k)
                  ssum(2) = ssum(2) + dpmid* saln(i,j,k,m)
                  ssum(3) = ssum(3) + dpnew* saln(i,j,k,n)
                endif !debug
!
                qdpmidn = 1.0/dpmidn
!
! ---           version that exactly conserves constant salinity
                smin   =   min( osaln(i,j,k), &
                                 saln(i,j,k,m), &
                                 saln(i,j,k,n) )
                dpsold = dpold*(osaln(i,j,k)  -smin)  !t-1
                dpsmid = dpmid*( saln(i,j,k,m)-smin)  !t
                dpsnew = dpnew*( saln(i,j,k,n)-smin)  !t+1
                q      = 0.5*ra2fac*(dpsold+dpsnew-2.0*dpsmid)
                saln(i,j,k,m)= smin + (dpsmid + q) * qdpmidn  !t with RA
!
                if     (ldebug_asselin .and. i.eq.itest &
                                       .and. j.eq.jtest) then
                  ssum(4) = ssum(4) + dpmidn*saln(i,j,k,m)
!                 write(6,*) k,ssum(4)*qonem
                  if     (k.eq.1) then
! ---               normalize by onem, layer 1 is often 1 m thick
                    s1(1) = ssum(1)/onem
                    s1(2) = ssum(2)/onem
                    s1(3) = ssum(3)/onem
                    s1(4) = ssum(4)/onem
                    write(lp,'(a,i9,1p4e16.8)') &
                     'asseln1:',nstep,s1(1),s1(2),s1(4),s1(3)
                  endif
                endif !debug
!
                if     (latemp) then
                  smin   =   min( otemp(i,j,k), &
                                   temp(i,j,k,m), &
                                   temp(i,j,k,n) )
                  dpsold = dpold*(otemp(i,j,k)  -smin)  !t-1
                  dpsmid = dpmid*( temp(i,j,k,m)-smin)  !t
                  dpsnew = dpnew*( temp(i,j,k,n)-smin)  !t+1
                  q      = 0.5*ra2fac*(dpsold+dpsnew-2.0*dpsmid)
                  temp(i,j,k,m) = smin + (dpsmid + q) * qdpmidn
! ---             update dependent thermodynamic variable
                  th3d(i,j,k,m) = sig(temp(i,j,k,m), &
                                      saln(i,j,k,m))-thbase
                elseif (lath3d) then
                  smin   =   min( oth3d(i,j,k), &
                                   th3d(i,j,k,m), &
                                   th3d(i,j,k,n) )
                  dpsold = dpold*(oth3d(i,j,k)  -smin)  !t-1
                  dpsmid = dpmid*( th3d(i,j,k,m)-smin)  !t
                  dpsnew = dpnew*( th3d(i,j,k,n)-smin)  !t+1
                  q      = 0.5*ra2fac*(dpsold+dpsnew-2.0*dpsmid)
                  th3d(i,j,k,m) = smin + (dpsmid + q) * qdpmidn
! ---             update dependent thermodynamic variable
                  temp(i,j,k,m) = tofsig(th3d(i,j,k,m)+thbase, &
                                         saln(i,j,k,m))
                else   ! exactly isopycnal layer
                  th3d(i,j,k,m) = theta(i,j,k)
! ---             update dependent thermodynamic variable
                  temp(i,j,k,m) = tofsig(th3d(i,j,k,m)+thbase, &
                                         saln(i,j,k,m))
                endif
                do ktr= 1,ntracr
! ---             non-negative version that exactly conserves constant tracers
                  smin   = min( otracer(i,j,k,  ktr), &
                                 tracer(i,j,k,m,ktr), &
                                 tracer(i,j,k,n,ktr) )
                  dpsold = dpold*(otracer(i,j,k,  ktr) - smin)  !>=0
                  dpsmid = dpmid*( tracer(i,j,k,m,ktr) - smin)  !>=0
                  dpsnew = dpnew*( tracer(i,j,k,n,ktr) - smin)  !>=0
                  q      = 0.5*ra2fac*(dpsold+dpsnew-2.0*dpsmid)
!diag               if (i.eq.itest.and.j.eq.jtest.and.ktr.eq.1) then
!diag                   write(lp,'(a,i3,1p3g15.6)') &
!diag                     'RA dpo:',k,dpold*qonem,dpmid*qonem,dpnew*qonem
!diag                   write(lp,'(a,i3,3f15.6)') &
!diag                     'RA tro:',k,otracer(i,j,k,  ktr), &
!diag                                  tracer(i,j,k,m,ktr), &
!diag                                  tracer(i,j,k,n,ktr)
!diag                   write(lp,'(a,i3,1p2g15.6)') &
!diag                     'RA dpm:',k,dpmid*qonem, &
!diag                                 dp(i,j,k,m)*qonem
!diag                   write(lp,'(a,i3,2f15.6)') &
!diag                     'RA trm:',k,tracer(i,j,k,m,ktr), &
!diag                                 smin + (dpsmid + q) * qdpmidn
!diag                   call flush(lp)
!diag               endif !debug
                  tracer(i,j,k,m,ktr) = smin + (dpsmid + q) * qdpmidn
                enddo !ktr
                if     (mxlmy) then
                  dpsold = dpold*oq2(i,j,k)
                  dpsmid = dpmid* q2(i,j,k,m)
                  dpsnew = dpnew* q2(i,j,k,n)
                  q      = 0.5*ra2fac*(dpsold+dpsnew-2.0*dpsmid)
                  q2(i,j,k,m) = (dpsmid + q) * qdpmidn
                  dpsold = dpold*oq2l(i,j,k)
                  dpsmid = dpmid* q2l(i,j,k,m)
                  dpsnew = dpnew* q2l(i,j,k,n)
                  q      = 0.5*ra2fac*(dpsold+dpsnew-2.0*dpsmid)
                  q2l(i,j,k,m) = (dpsmid + q) * qdpmidn
                endif !mxlmy
              endif !effectively non-zero 
!
              if     (ldebug_asselin .and. k.eq.kk .and. &
                          i.eq.itest .and. j.eq.jtest   ) then
! ---           normalize by initial depth
                ssum(1) = ssum(1)/pbot(i,j)
                ssum(2) = ssum(2)/pbot(i,j)
                ssum(3) = ssum(3)/pbot(i,j)
                ssum(4) = ssum(4)/pbot(i,j)
                write(lp,'(a,i9,1p4e16.8)') &
                 'asselin:',nstep,ssum(1),ssum(2),ssum(4),ssum(3)
              endif !debug
!
            endif !ip
          enddo !i
        enddo !k
      enddo !j
!$OMP END PARALLEL DO
!
      if     (lpipe .and. lpipe_asselin) then
! ---   compare two model runs.
        do k= 1,kk
          write (text,'(a9,i3)') 'otemp  k=',k
          call pipe_compare_sym1(otemp(1-nbdy,1-nbdy,k),  ip,text)
          write (text,'(a9,i3)') 'temp.m k=',k
          call pipe_compare_sym1( temp(1-nbdy,1-nbdy,k,m),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)') 'osaln  k=',k
          call pipe_compare_sym1(osaln(1-nbdy,1-nbdy,k),  ip,text)
          write (text,'(a9,i3)') 'saln.m k=',k
          call pipe_compare_sym1( saln(1-nbdy,1-nbdy,k,m),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)') 'oth3d  k=',k
          call pipe_compare_sym1(oth3d(1-nbdy,1-nbdy,k),  ip,text)
          write (text,'(a9,i3)') 'th3d.m k=',k
          call pipe_compare_sym1( th3d(1-nbdy,1-nbdy,k,m),ip,text)
          write (text,'(a9,i3)') 'th3d.n k=',k
          call pipe_compare_sym1( th3d(1-nbdy,1-nbdy,k,n),ip,text)
        enddo !k
      endif
!
!     call pipe_comparall(m,n, 'asseln, step')
!
      return
      end subroutine asselin_filter

      end module mod_asselin
!
!
!> Revision history:
!>
!> Aug. 2018 - 1st version
!> Feb. 2019 - replaced onetai by 1.0
!> Sep. 2019 - added oneta0