c-----------------------------------------------------------------------
      subroutine shapu_diff(dt,npt,nz,u,v,h,fu,fv,lxxk,lyxk,isyk,isk,
     *           ifxk,ifpxk,ifyk,tp)
c-----------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_para.h'
      dimension u(npt,1),v(npt,1),h(npt,1),tp(1),fu(npt,1),fv(npt,1)
      dimension lxxk(MXBDY,nz),lyxk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *          ifxk(9*MAXSID,nz), ifpxk(5*MAXSID,nz), ifyk(9*MAXSID,nz)
      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +         ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      logical NOSLIP
      common /new_shap/nordu,nshapu,mshapu,dshapu, nordh,nshaph,mshaph,dshaph
     *               , nordp,nshapp,mshapp,dshapp
     *               , nordd,nshapd,mshapd,dshapd,dshapv
      common/param0/iyear,iday,isec,delt,ncyc,mbc,nonlin,label(20),
     +              itherm,mlc,limp
      common /new_filt/  MAXFO, nbx, nby, nfx, nfpx, nfy

      external shap_1drn, shap_1dro

      if (mshapu .eq. 0) return

      if (nshapu.ne.0) then
         NOSLIP = (mbc .eq. 1)
         do k = 1, nz
            npk = nptk(k) 
            nbx = nbxk(k) 
            nby = nbyk(k) 
            nfx = nfxk(k) 
            nfpx = nfpxk(k) 
            nfy = nfyk(k) 
            if     (mshapu .eq. 1) then
c........use NEW reduce order filter near BC & in NARROW passages: 
               call dshap_2d (dt,nordu,.true.,.true.,NOSLIP,npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),u(1,k),h(1,k),fu(1,k),tp,shap_1drn,dshapu,dshapv)
               call dshap_2d (dt,nordu,.true.,NOSLIP,.true., npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),v(1,k),h(1,k),fv(1,k),tp,shap_1drn,dshapv,dshapu)
            elseif (mshapu .eq. 2) then
c........use NEW reduce order filter near BC only:
               call dshap_2d (dt,nordu,.false.,.true.,NOSLIP,npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),u(1,k),h(1,k),fu(1,k),tp,shap_1drn,dshapu,dshapv)
               call dshap_2d (dt,nordu,.false.,NOSLIP,.true., npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),v(1,k),h(1,k),fv(1,k),tp,shap_1drn,dshapv,dshapu)
            elseif (mshapu .eq. 3) then
c........use OLD reduce order filter near BC & in NARROW passages: 
               call dshap_2d (dt,nordu,.true.,.true.,NOSLIP,npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),u(1,k),h(1,k),fu(1,k),tp,shap_1dro,dshapu,dshapv)
               call dshap_2d (dt,nordu,.true.,NOSLIP,.true., npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),v(1,k),h(1,k),fv(1,k),tp,shap_1dro,dshapv,dshapu)
            endif
         enddo
      endif
            
      return
      end


c--------------------------------------------------------------------------
      subroutine dshap_2d (dt,nord,REDUC,SETX,SETY,npt,lxx,lyx,isy,isk,
     *                     ifx,ifpx,ify, aa,h,ff,abc,filter,xcoef,ycoef)
c--------------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension lxx(1),lyx(1),isy(1),isk(1), aa(1), bb(1), cc(1), abc(npt,1)
      dimension ifx(3,1), ifpx(5,1), ify(3,1) , ff(1), h(1)
      logical SETX, SETY, ZER, REDUC 
      pointer (pbb, bb), (pcc, cc)
      external filter

      common /shap_c25/ cs25(10), ZER, s_coef
      common /new_filt/  MAXFO, nbx, nby, nfx, nfpx, nfy

      s_coef = xcoef

c sizeof(abc) = 3*npt

      pbb = loc(abc(1,1))
      pcc = loc(abc(1,2))


c.......apply Shapiro filter in X direction
      do i = 1, npt
         bb(i) = aa(isk(i))
      enddo

      ZER = SETX   ! used only if filter = shap_1do
      nshap = MAXFO

      do k = 1, nfx
         if (REDUC) nshap = min(nord, ifx(3,k))
         call filter (nshap, ifx(2,k), bb(ifx(1,k)), cc)
      enddo

      do k = 1, nfpx
         if (REDUC) nshap = min(nord, ifpx(5,k))
         if (ifpx(4,k) .eq. 0) then
            call shap_1dp0 (nshap, ifpx(2,k), bb(ifpx(1,k)), cc)
         else
            call shap_1dper (nshap, ifpx(1,k), bb, cc, filter)
         endif
      enddo

      do i = 1, npt
         cc(isk(i)) = bb(i)
      enddo
      
      if (SETX) then
         do i = 1, nbx
            cc(lxx(i)) = 0.
         enddo
      endif
      if (SETY) then
         do i = 1, nby
            cc(lyx(i)) = 0.
         enddo
      endif

      do j = 1, npt
         i = isk(j)
         ff(i) = ff(i) - h(i)*cc(i)/dt
      enddo
      
c.......apply Shapiro filter in Y direction
      s_coef = ycoef

      do i = 1, npt
         bb(i) = aa(isy(i))
      enddo

      ZER = SETY 
      do k = 1, nfy
         if (REDUC) nshap = min(nord, ify(3,k))
         call filter (nshap, ify(2,k), bb(ify(1,k)), cc)
      enddo

      do i = 1, npt
         cc(isy(i)) = bb(i)
      enddo

      if (SETX) then
         do i = 1, nbx
            cc(lxx(i)) = 0.
         enddo
      endif
      if (SETY) then
         do i = 1, nby
            cc(lyx(i)) = 0.
         enddo
      endif

      do j = 1, npt
         i = isk(j)
         ff(i) = ff(i) - h(i)*cc(i)/dt
      enddo

      return
      end
c-----------------------------------------------------------------------
      subroutine shapuh_diff(dt,npt,nsig,nz,u,v,h,fu,fv,lxxk,lyxk,isyk,
     *           isk,ifxk,ifpxk,ifyk,tp)
c-----------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_para.h'
      dimension u(npt,1),v(npt,1),h(npt,1),tp(1),fu(npt,1),fv(npt,1)
      dimension lxxk(MXBDY,nz),lyxk(MAXNB,nz),isyk(npt,nz),isk(npt,nz),
     *          ifxk(9*MAXSID,nz), ifpxk(5*MAXSID,nz), ifyk(9*MAXSID,nz)
      common/gridk/nptk(MAXNZ),nbxk(MAXNZ),nbyk(MAXNZ),ncsk(MAXNZ)
     +         ,npbck(MAXNZ),nlok(MAXNZ),nfxk(MAXNZ),nfpxk(MAXNZ),nfyk(MAXNZ)
      logical NOSLIP
      common /new_shap/nordu,nshapu,mshapu,dshapu, nordh,nshaph,mshaph,dshaph
     *               , nordp,nshapp,mshapp,dshapp
     *               , nordd,nshapd,mshapd,dshapd,dshapv
      common/param0/iyear,iday,isec,delt,ncyc,mbc,nonlin,label(20),
     +              itherm,mlc,limp
      common /new_filt/  MAXFO, nbx, nby, nfx, nfpx, nfy

      external shap_1drn, shap_1dro

      if (mshapu .eq. 0) return

      if (nshapu.ne.0) then
         NOSLIP = (mbc .eq. 1)
         do k = 1, nsig
            npk = nptk(k) 
            nbx = nbxk(k) 
            nby = nbyk(k) 
            nfx = nfxk(k) 
            nfpx = nfpxk(k) 
            nfy = nfyk(k) 
            if     (mshapu .eq. 1) then
c........use NEW reduce order filter near BC & in NARROW passages: 
               call dshap_2dh (dt,nordu,.true.,.true.,NOSLIP,npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),u(1,k),h(1,k),fu(1,k),tp,shap_1drn,dshapu,dshapv)
               call dshap_2dh (dt,nordu,.true.,NOSLIP,.true., npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),v(1,k),h(1,k),fv(1,k),tp,shap_1drn,dshapv,dshapu)
            elseif (mshapu .eq. 2) then
c........use NEW reduce order filter near BC only:
               call dshap_2dh (dt,nordu,.false.,.true.,NOSLIP,npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),u(1,k),h(1,k),fu(1,k),tp,shap_1drn,dshapu,dshapv)
               call dshap_2dh (dt,nordu,.false.,NOSLIP,.true., npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),v(1,k),h(1,k),fv(1,k),tp,shap_1drn,dshapv,dshapu)
            elseif (mshapu .eq. 3) then
c........use OLD reduce order filter near BC & in NARROW passages: 
               call dshap_2dh (dt,nordu,.true.,.true.,NOSLIP,npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),u(1,k),h(1,k),fu(1,k),tp,shap_1dro,dshapu,dshapv)
               call dshap_2dh (dt,nordu,.true.,NOSLIP,.true., npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),v(1,k),h(1,k),fv(1,k),tp,shap_1dro,dshapv,dshapu)
            endif
         enddo
         do k = nsig+1, nz
            npk = nptk(k) 
            nbx = nbxk(k) 
            nby = nbyk(k) 
            nfx = nfxk(k) 
            nfpx = nfpxk(k) 
            nfy = nfyk(k) 
            if     (mshapu .eq. 1) then
c........use NEW reduce order filter near BC & in NARROW passages: 
               call dshap_2d (dt,nordu,.true.,.true.,NOSLIP,npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),u(1,k),h(1,k),fu(1,k),tp,shap_1drn,dshapu,dshapv)
               call dshap_2d (dt,nordu,.true.,NOSLIP,.true., npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),v(1,k),h(1,k),fv(1,k),tp,shap_1drn,dshapv,dshapu)
            elseif (mshapu .eq. 2) then
c........use NEW reduce order filter near BC only:
               call dshap_2d (dt,nordu,.false.,.true.,NOSLIP,npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),u(1,k),h(1,k),fu(1,k),tp,shap_1drn,dshapu,dshapv)
               call dshap_2d (dt,nordu,.false.,NOSLIP,.true., npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),v(1,k),h(1,k),fv(1,k),tp,shap_1drn,dshapv,dshapu)
            elseif (mshapu .eq. 3) then
c........use OLD reduce order filter near BC & in NARROW passages: 
               call dshap_2d (dt,nordu,.true.,.true.,NOSLIP,npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),u(1,k),h(1,k),fu(1,k),tp,shap_1dro,dshapu,dshapv)
               call dshap_2d (dt,nordu,.true.,NOSLIP,.true., npk, 
     *           lxxk(1,k),lyxk(1,k),isyk(1,k),isk(1,k),ifxk(1,k),ifpxk(1,k)
     *           ,ifyk(1,k),v(1,k),h(1,k),fv(1,k),tp,shap_1dro,dshapv,dshapu)
            endif
         enddo
      endif
            
      return
      end


c--------------------------------------------------------------------------
      subroutine dshap_2dh (dt,nord,REDUC,SETX,SETY,npt,lxx,lyx,isy,isk,
     *                     ifx,ifpx,ify, aa,h,ff,abc,filter,scoef)
c--------------------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      dimension lxx(1),lyx(1),isy(1),isk(1), aa(1), bb(1), cc(1), abc(npt,1)
      dimension ifx(3,1), ifpx(5,1), ify(3,1) , ff(1), h(1)
      logical SETX, SETY, ZER, REDUC 
      pointer (pbb, bb), (pcc, cc)
      external filter

      common /shap_c25/ cs25(10), ZER, s_coef
      common /new_filt/  MAXFO, nbx, nby, nfx, nfpx, nfy

      s_coef = scoef

c sizeof(abc) = 3*npt

      pbb = loc(abc(1,1))
      pcc = loc(abc(1,2))

c.......apply Shapiro filter in X direction
      do j = 1, npt
         i = isk(j)
         bb(j) = h(i)*aa(i)
      enddo

      ZER = SETX   ! used only if filter = shap_1do
      nshap = MAXFO
      do k = 1, nfx
         if (REDUC) nshap = min(nord, ifx(3,k))
         call filter (nshap, ifx(2,k), bb(ifx(1,k)), cc)
      enddo

      do k = 1, nfpx
         if (REDUC) nshap = min(nord, ifpx(5,k))
         if (ifpx(4,k) .eq. 0) then
            call shap_1dp0 (nshap, ifpx(2,k), bb(ifpx(1,k)), cc)
         else
            call shap_1dper (nshap, ifpx(1,k), bb, cc, filter)
         endif
      enddo

      do i = 1, npt
         cc(isk(i)) = bb(i)
      enddo
      
      if (SETX) then
         do i = 1, nbx
            cc(lxx(i)) = 0.
         enddo
      endif
      if (SETY) then
         do i = 1, nby
            cc(lyx(i)) = 0.
         enddo
      endif

      do j = 1, npt
         i = isk(j)
         ff(i) = ff(i) - cc(i)/dt
      enddo
      
c.......apply Shapiro filter in Y direction
      do j = 1, npt
         i = isy(j)
         bb(j) = h(i)*aa(i)
      enddo

      ZER = SETY 
      do k = 1, nfy
         if (REDUC) nshap = min(nord, ify(3,k))
         call filter (nshap, ify(2,k), bb(ify(1,k)), cc)
      enddo

      do i = 1, npt
         cc(isy(i)) = bb(i)
      enddo

      if (SETX) then
         do i = 1, nbx
            cc(lxx(i)) = 0.
         enddo
      endif
      if (SETY) then
         do i = 1, nby
            cc(lyx(i)) = 0.
         enddo
      endif

      do j = 1, npt
         i = isk(j)
         ff(i) = ff(i) - cc(i)/dt
      enddo

      return
      end
