C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/shap_filt/shap_filt_uv_s4.F,v 1.2 2001/05/29 14:01:40 adcroft Exp $
C $Name:  $

#include "SHAP_FILT_OPTIONS.h"

      SUBROUTINE SHAP_FILT_UV_S4( 
     U           uFld, vFld,
     I           myTime, myThid )
C     /==========================================================\
C     | S/R SHAP_FILT_UV                                         |
C     | Applies Shapiro filter to tracer field over one XY slice |
C     | of one tile at a time.                                   |
C     \==========================================================/
      IMPLICIT NONE

C     == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_SHAP_FILT
#include "SHAP_FILT.h"
#include "SHAP_FILT_UV.h"
#endif

C     == Routine arguments
      _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL     myTime
      INTEGER myThid

#ifdef ALLOW_SHAP_FILT

C     == Local variables ==
      INTEGER bi,bj,K,I,J,N
      _RL tmpGrdU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL tmpGrdV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL maskZj,maskZp

      IF (nShapUV.gt.0) THEN

        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO K=1,Nr
           DO J=1,sNy
            DO I=1,sNx+1
             tmpFldU(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)
     &                *_maskW(i,j,k,bi,bj)
            ENDDO
           ENDDO
           DO J=1,sNy+1
            DO I=1,sNx
             tmpFldV(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)
     &                *_maskS(i,j,k,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO


C      d_xx^n tmpFld 

       DO N=1,nShapUV

        CALL EXCH_UV_XYZ_RL(tmpFldU,tmpFldV,.TRUE.,myThid)


        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO K=1,Nr

C          Uxx
           DO J=1,sNy
            DO I=1,sNx+1
             tmpGrdU(i,j) = -0.25*(
     &          tmpFldU(i-1,j,k,bi,bj) + tmpFldU(i+1,j,k,bi,bj)
     &             - 2.*tmpFldU(i,j,k,bi,bj)
     &            )*_maskW(i,j,k,bi,bj)
            ENDDO
           ENDDO

           DO J=1,sNy
            DO I=1,sNx+1
             tmpFldU(i,j,k,bi,bj) = tmpGrdU(i,j)
            ENDDO
           ENDDO

C          Vyy
           DO J=1,sNy+1
            DO I=1,sNx
             tmpGrdV(i,j) = -0.25*(
     &          tmpFldV(i,j-1,k,bi,bj) + tmpFldV(i,j+1,k,bi,bj)
     &             - 2.*tmpFldV(i,j,k,bi,bj)
     &            )*_maskS(i,j,k,bi,bj)
            ENDDO
           ENDDO

           DO J=1,sNy+1
            DO I=1,sNx
             tmpFldV(i,j,k,bi,bj) = tmpGrdV(i,j)
            ENDDO
           ENDDO

          ENDDO
         ENDDO
        ENDDO

       ENDDO

C      F <-  [1-d_xx^n]F
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO K=1,Nr
          DO J=1,sNy
           DO I=1,sNx+1
            tmpFldU(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)-tmpFldU(i,j,k,bi,bj)
            uFld(i,j,k,bi,bj)=tmpFldU(i,j,k,bi,bj)
           ENDDO
          ENDDO
          DO J=1,sNy+1
           DO I=1,sNx
            tmpFldV(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)-tmpFldV(i,j,k,bi,bj)
            vFld(i,j,k,bi,bj)=tmpFldV(i,j,k,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO


C      d_yy^n tmpFld 

       DO N=1,nShapUV

        CALL EXCH_UV_XYZ_RL(tmpFldU,tmpFldV,.TRUE.,myThid)

        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO K=1,Nr

C          Uyy
           DO J=1,sNy
            DO I=1,sNx+1
             maskZj=_maskS(i-1, j ,k,bi,bj)
     &             *_maskS( i , j ,k,bi,bj)
             maskZp=_maskS(i-1,j+1,k,bi,bj)
     &             *_maskS( i ,j+1,k,bi,bj)
             tmpGrdU(i,j) = -0.25*(
     &        (tmpFldU(i,j+1,k,bi,bj)-tmpFldU(i, j ,k,bi,bj))*maskZp
     &       -(tmpFldU(i, j ,k,bi,bj)-tmpFldU(i,j-1,k,bi,bj))*maskZj
#ifdef NO_SLIP_SHAP
     &       -2.*(2.-maskZj-maskZp)*tmpFldU(i,j,k,bi,bj)
#endif
     &             )*_maskW(i,j,k,bi,bj)
            ENDDO
           ENDDO

           IF (useCubedSphereExchange) THEN
            J=1
            DO I=1,sNx+1,sNx
             maskZj=maskS(i-1, j ,k,bi,bj)*maskS( i , j ,k,bi,bj)
             maskZp=maskS(i-1,j+1,k,bi,bj)*maskS( i ,j+1,k,bi,bj)
             tmpGrdU(i,j) = -0.25*(
     &        (tmpFldU(i,j+1,k,bi,bj)-tmpFldU(i, j ,k,bi,bj))*maskZp
     &       -(tmpFldU(i, j ,k,bi,bj)-0*tmpFldU(i,j-1,k,bi,bj))*maskZj
#ifdef NO_SLIP_SHAP
     &       -2.*(2.-maskZj-maskZp)*tmpFldU(i,j,k,bi,bj)
#endif
     &             )*_maskW(i,j,k,bi,bj)
            ENDDO
            J=sNy
            DO I=1,sNx+1,sNx
             maskZj=maskS(i-1, j ,k,bi,bj)*maskS( i , j ,k,bi,bj)
             maskZp=maskS(i-1,j+1,k,bi,bj)*maskS( i ,j+1,k,bi,bj)
             tmpGrdU(i,j) = -0.25*(
     &        (0*tmpFldU(i,j+1,k,bi,bj)-tmpFldU(i, j ,k,bi,bj))*maskZp
     &       -(tmpFldU(i, j ,k,bi,bj)-tmpFldU(i,j-1,k,bi,bj))*maskZj
#ifdef NO_SLIP_SHAP
     &       -2.*(2.-maskZj-maskZp)*tmpFldU(i,j,k,bi,bj)
#endif
     &             )*_maskW(i,j,k,bi,bj)
            ENDDO
           ENDIF

           DO J=1,sNy
            DO I=1,sNx+1
             tmpFldU(i,j,k,bi,bj) = tmpGrdU(i,j)
            ENDDO
           ENDDO

C          Vxx
           DO J=1,sNy+1
            DO I=1,sNx
             maskZj=_maskW( i ,j-1,k,bi,bj)
     &             *_maskW( i , j ,k,bi,bj)
             maskZp=_maskW(i+1,j-1,k,bi,bj)
     &             *_maskW(i+1, j ,k,bi,bj)
             tmpGrdV(i,j) = -0.25*(
     &        (tmpFldV(i+1,j,k,bi,bj)-tmpFldV( i ,j,k,bi,bj))*maskZp
     &       -(tmpFldV( i ,j,k,bi,bj)-tmpFldV(i-1,j,k,bi,bj))*maskZj
#ifdef NO_SLIP_SHAP
     &       -2.*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
#endif
     &             )*_maskS(i,j,k,bi,bj)
            ENDDO
           ENDDO

           IF (useCubedSphereExchange) THEN
            DO J=1,sNy+1,sNy
            I=1
             maskZj=maskW( i ,j-1,k,bi,bj)*maskW( i , j ,k,bi,bj)
             maskZp=maskW(i+1,j-1,k,bi,bj)*maskW(i+1, j ,k,bi,bj)
             tmpGrdV(i,j) = -0.25*(
     &        (tmpFldV(i+1,j,k,bi,bj)-tmpFldV( i ,j,k,bi,bj))*maskZp
     &       -(tmpFldV( i ,j,k,bi,bj)-0*tmpFldV(i-1,j,k,bi,bj))*maskZj
#ifdef NO_SLIP_SHAP
     &       -2.*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
#endif
     &             )*_maskS(i,j,k,bi,bj)
            ENDDO
            DO J=1,sNy+1,sNy
            I=sNx
             maskZj=maskW( i ,j-1,k,bi,bj)*maskW( i , j ,k,bi,bj)
             maskZp=maskW(i+1,j-1,k,bi,bj)*maskW(i+1, j ,k,bi,bj)
             tmpGrdV(i,j) = -0.25*(
     &        (0*tmpFldV(i+1,j,k,bi,bj)-tmpFldV( i ,j,k,bi,bj))*maskZp
     &       -(tmpFldV( i ,j,k,bi,bj)-tmpFldV(i-1,j,k,bi,bj))*maskZj
#ifdef NO_SLIP_SHAP
     &       -2.*(2.-maskZj-maskZp)*tmpFldV(i,j,k,bi,bj)
#endif
     &             )*_maskS(i,j,k,bi,bj)
            ENDDO
           ENDIF

           DO J=1,sNy+1
            DO I=1,sNx
             tmpFldV(i,j,k,bi,bj) = tmpGrdV(i,j)
            ENDDO
           ENDDO

          ENDDO
         ENDDO
        ENDDO

       ENDDO

C      F <-  [1-d_yy^n]F
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO K=1,Nr
          DO J=1,sNy
           DO I=1,sNx+1
            uFld(i,j,k,bi,bj)=uFld(i,j,k,bi,bj)-tmpFldU(i,j,k,bi,bj)
           ENDDO
          ENDDO
          DO J=1,sNy+1
           DO I=1,sNx
            vFld(i,j,k,bi,bj)=vFld(i,j,k,bi,bj)-tmpFldV(i,j,k,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO

      ENDIF
#endif /* ALLOW_SHAP_FILT */

      RETURN
      END