pkg/shap_filt/shap_filt_tracer_s2.F

C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.8 2005/09/27 22:11:06 jmc Exp $
C $Name:  $

#include "SHAP_FILT_OPTIONS.h"

CBOP
C     !ROUTINE: SHAP_FILT_TRACER_S2
C     !INTERFACE:
      SUBROUTINE SHAP_FILT_TRACER_S2(
     U           field, tmpFld,
     I           nShapTr, kSize, myTime, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R SHAP_FILT_TRACER_S2
C     | o Applies Shapiro filter to 2D field (cell center).
C     | o use filtering function "S2" = [1 - (d_xx+d_yy)^n]
C     | o Options for computational filter (no grid spacing)
C     |   or physical space filter (with grid spacing) or both.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SHAP_FILT.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments
C     field :: cell-centered 2D field on which filter applies
C     tmpFld :: working temporary array
C     nShapTr :: (total) power of the filter for this tracer
C     kSize :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
C     myTime :: Current time in simulation
C     myThid :: Thread number for this instance of SHAP_FILT_TRACER_S2
      INTEGER nShapTr, kSize
      _RL field(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
      _RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
      _RL     myTime
      INTEGER myThid

#ifdef ALLOW_SHAP_FILT

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER nShapComput
      INTEGER bi,bj,k,i,j,n
      _RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL tmpFdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL tmpFdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

      IF (nShapTr.GT.0) THEN
C-------
C  Apply computational filter ^(nShap-nShapPhys) without grid factor
C  then apply Physical filter ^nShapPhys  with grid factors
C-------
        nShapComput = nShapTr - nShapTrPhys

        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO k=1,kSize
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO


C      ( d_xx +d_yy )^n tmpFld

       DO n=1,nShapTr

        IF ( MOD(n,2).EQ.1 .OR. Shap_alwaysExchTr ) THEN
         IF (kSize.EQ.Nr) THEN
          _EXCH_XYZ_R8( tmpFld, myThid )
         ELSEIF (kSize.EQ.1) THEN
          _EXCH_XY_R8( tmpFld, myThid )
         ELSE
          STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
         ENDIF
        ENDIF

        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO k=1,kSize

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--        Calculate gradient in X direction:
           IF ( .NOT.Shap_alwaysExchTr
     &          .AND. useCubedSphereExchange ) THEN
C          to compute d/dx(tmpFld), fill corners with appropriate values:
             CALL FILL_CS_CORNER_TR_RL( .TRUE., 
     &                                 tmpFld(1-Olx,1-Oly,k,bi,bj), 
     &                                 bi,bj, myThid )
           ENDIF
           IF ( n.LE.nShapComput ) THEN
C-         Computational space: del_i
             DO j=0,sNy+1
              DO i=0,sNx+2
               tmpFdx(i,j) = 
     &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
     &                     *_maskW(i,j,k,bi,bj)
              ENDDO
             ENDDO
           ELSE
C-         Physical space: grad_x
             DO j=0,sNy+1
              DO i=0,sNx+2
               tmpFdx(i,j) = 
     &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
     &                     *_hFacW(i,j,k,bi,bj)
     &                     *dyG(i,j,bi,bj)*recip_dxC(i,j,bi,bj)
              ENDDO
             ENDDO
           ENDIF

C--        Calculate gradient in Y direction:
           IF ( .NOT.Shap_alwaysExchTr
     &          .AND. useCubedSphereExchange ) THEN
C          to compute d/dy(tmpFld), fill corners with appropriate values:
             CALL FILL_CS_CORNER_TR_RL(.FALSE., 
     &                                 tmpFld(1-Olx,1-Oly,k,bi,bj), 
     &                                 bi,bj, myThid )
           ENDIF
           IF ( n.LE.nShapComput ) THEN
C-         Computational space: del_j
             DO j=0,sNy+2
              DO i=0,sNx+1
               tmpFdy(i,j) = 
     &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
     &                     *_maskS(i,j,k,bi,bj)
              ENDDO
             ENDDO
           ELSE
C-         Physical space: grad_y
             DO j=0,sNy+2
              DO i=0,sNx+1
               tmpFdy(i,j) = 
     &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
     &                     *_hFacS(i,j,k,bi,bj)
     &                     *dxG(i,j,bi,bj)*recip_dyC(i,j,bi,bj)
              ENDDO
             ENDDO
           ENDIF

C--        Calculate (d_xx + d_yy) tmpFld :
           DO j=0,sNy+1
             DO i=0,sNx+1
               tmpGrd(i,j) = ( tmpFdx(i+1,j) - tmpFdx(i,j) )
     &                     + ( tmpFdy(i,j+1) - tmpFdy(i,j) )
             ENDDO
           ENDDO

C--        Computational space Filter 
           IF ( n.LE.nShapComput ) THEN
             DO j=0,sNy+1
              DO i=0,sNx+1
               tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
              ENDDO
             ENDDO
C--        Physical space Filter 
           ELSEIF (Shap_TrLength.LE.0.) THEN
             DO j=0,sNy+1
              DO i=0,sNx+1
               tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
     &             *recip_hFacC(i,j,k,bi,bj)
              ENDDO
             ENDDO
           ELSE
             DO j=0,sNy+1
              DO i=0,sNx+1
               tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
     &             *recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
     &             *Shap_TrLength*Shap_TrLength
              ENDDO
             ENDDO
           ENDIF
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C      end k,bi,bj loop:
          ENDDO
         ENDDO
        ENDDO
C      end loop n=1,nShapTr
       ENDDO

C      F <-  [1 - (d_xx+d_yy)^n *deltaT/tau].F
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO k=1,kSize
          DO j=1,sNy
           DO i=1,sNx
            field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
     &            -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
            tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO

c       IF (kSize.EQ.Nr) THEN
c         _EXCH_XYZ_R8( field, myThid )
c       ELSEIF (kSize.EQ.1) THEN
c         _EXCH_XY_R8( field, myThid )
c       ELSE
c         STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
c       ENDIF

      ENDIF
#endif /* ALLOW_SHAP_FILT */

      RETURN
      END
1	jmc	1.9	C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.8 2005/09/27 22:11:06 jmc Exp $
2	jmc	1.3	C $Name: $
3	adcroft	1.2
4			#include "SHAP_FILT_OPTIONS.h"
5	jmc	1.8
6	jmc	1.4	CBOP
7			C !ROUTINE: SHAP_FILT_TRACER_S2
8			C !INTERFACE:
9			SUBROUTINE SHAP_FILT_TRACER_S2(
10			U field, tmpFld,
11	jmc	1.5	I nShapTr, kSize, myTime, myThid )
12	jmc	1.4	C !DESCRIPTION: \bv
13			C ==========================================================
14			C \| S/R SHAP_FILT_TRACER_S2
15			C \| o Applies Shapiro filter to 2D field (cell center).
16			C \| o use filtering function "S2" = [1 - (d_xx+d_yy)^n]
17			C \| o Options for computational filter (no grid spacing)
18			C \| or physical space filter (with grid spacing) or both.
19			C ==========================================================
20			C \ev
21	jmc	1.8
22	jmc	1.4	C !USES:
23	adcroft	1.2	IMPLICIT NONE
24	jmc	1.8
25	adcroft	1.2	C == Global variables ===
26			#include "SIZE.h"
27			#include "EEPARAMS.h"
28			#include "PARAMS.h"
29			#include "GRID.h"
30			#include "SHAP_FILT.h"
31
32	jmc	1.4	C !INPUT/OUTPUT PARAMETERS:
33	adcroft	1.2	C == Routine arguments
34	jmc	1.4	C field :: cell-centered 2D field on which filter applies
35			C tmpFld :: working temporary array
36	jmc	1.5	C nShapTr :: (total) power of the filter for this tracer
37	jmc	1.4	C kSize :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
38			C myTime :: Current time in simulation
39			C myThid :: Thread number for this instance of SHAP_FILT_TRACER_S2
40	jmc	1.5	INTEGER nShapTr, kSize
41	jmc	1.4	_RL field(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
42			_RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
43	adcroft	1.2	_RL myTime
44			INTEGER myThid
45	jmc	1.8
46	adcroft	1.2	#ifdef ALLOW_SHAP_FILT
47
48	jmc	1.4	C !LOCAL VARIABLES:
49	adcroft	1.2	C == Local variables ==
50	jmc	1.9	INTEGER nShapComput
51			INTEGER bi,bj,k,i,j,n
52	adcroft	1.2	_RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53	jmc	1.9	_RL tmpFdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54			_RL tmpFdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55	jmc	1.4	CEOP
56	adcroft	1.2
57	jmc	1.8	IF (nShapTr.GT.0) THEN
58	jmc	1.3	C-------
59			C Apply computational filter ^(nShap-nShapPhys) without grid factor
60	jmc	1.9	C then apply Physical filter ^nShapPhys with grid factors
61	jmc	1.3	C-------
62	jmc	1.9	nShapComput = nShapTr - nShapTrPhys
63	adcroft	1.2
64			DO bj=myByLo(myThid),myByHi(myThid)
65			DO bi=myBxLo(myThid),myBxHi(myThid)
66	jmc	1.9	DO k=1,kSize
67			DO j=1-Oly,sNy+Oly
68			DO i=1-Olx,sNx+Olx
69	adcroft	1.2	tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
70			ENDDO
71			ENDDO
72			ENDDO
73			ENDDO
74			ENDDO
75
76
77	jmc	1.9	C ( d_xx +d_yy )^n tmpFld
78	adcroft	1.2
79	jmc	1.9	DO n=1,nShapTr
80	adcroft	1.2
81	jmc	1.9	IF ( MOD(n,2).EQ.1 .OR. Shap_alwaysExchTr ) THEN
82			IF (kSize.EQ.Nr) THEN
83	jmc	1.4	_EXCH_XYZ_R8( tmpFld, myThid )
84	jmc	1.9	ELSEIF (kSize.EQ.1) THEN
85	jmc	1.4	_EXCH_XY_R8( tmpFld, myThid )
86	jmc	1.9	ELSE
87			STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
88			ENDIF
89	jmc	1.4	ENDIF
90	jmc	1.8
91	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
92			DO bi=myBxLo(myThid),myBxHi(myThid)
93	jmc	1.9	DO k=1,kSize
94	jmc	1.8
95	jmc	1.9	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
96	adcroft	1.2
97	jmc	1.9	C-- Calculate gradient in X direction:
98			IF ( .NOT.Shap_alwaysExchTr
99			& .AND. useCubedSphereExchange ) THEN
100			C to compute d/dx(tmpFld), fill corners with appropriate values:
101			CALL FILL_CS_CORNER_TR_RL( .TRUE.,
102			& tmpFld(1-Olx,1-Oly,k,bi,bj),
103			& bi,bj, myThid )
104			ENDIF
105			IF ( n.LE.nShapComput ) THEN
106			C- Computational space: del_i
107			DO j=0,sNy+1
108			DO i=0,sNx+2
109			tmpFdx(i,j) =
110			& ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
111			& *_maskW(i,j,k,bi,bj)
112			ENDDO
113			ENDDO
114			ELSE
115			C- Physical space: grad_x
116			DO j=0,sNy+1
117			DO i=0,sNx+2
118			tmpFdx(i,j) =
119			& ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
120			& *_hFacW(i,j,k,bi,bj)
121			& dyG(i,j,bi,bj)recip_dxC(i,j,bi,bj)
122			ENDDO
123			ENDDO
124			ENDIF
125	adcroft	1.2
126	jmc	1.9	C-- Calculate gradient in Y direction:
127			IF ( .NOT.Shap_alwaysExchTr
128			& .AND. useCubedSphereExchange ) THEN
129			C to compute d/dy(tmpFld), fill corners with appropriate values:
130			CALL FILL_CS_CORNER_TR_RL(.FALSE.,
131			& tmpFld(1-Olx,1-Oly,k,bi,bj),
132			& bi,bj, myThid )
133			ENDIF
134			IF ( n.LE.nShapComput ) THEN
135			C- Computational space: del_j
136			DO j=0,sNy+2
137			DO i=0,sNx+1
138			tmpFdy(i,j) =
139			& ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
140			& *_maskS(i,j,k,bi,bj)
141			ENDDO
142			ENDDO
143			ELSE
144			C- Physical space: grad_y
145			DO j=0,sNy+2
146			DO i=0,sNx+1
147			tmpFdy(i,j) =
148			& ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
149			& *_hFacS(i,j,k,bi,bj)
150			& dxG(i,j,bi,bj)recip_dyC(i,j,bi,bj)
151			ENDDO
152			ENDDO
153			ENDIF
154	jmc	1.3
155	jmc	1.9	C-- Calculate (d_xx + d_yy) tmpFld :
156			DO j=0,sNy+1
157			DO i=0,sNx+1
158			tmpGrd(i,j) = ( tmpFdx(i+1,j) - tmpFdx(i,j) )
159			& + ( tmpFdy(i,j+1) - tmpFdy(i,j) )
160			ENDDO
161	jmc	1.3	ENDDO
162
163	jmc	1.9	C-- Computational space Filter
164			IF ( n.LE.nShapComput ) THEN
165			DO j=0,sNy+1
166			DO i=0,sNx+1
167			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
168			ENDDO
169			ENDDO
170			C-- Physical space Filter
171			ELSEIF (Shap_TrLength.LE.0.) THEN
172			DO j=0,sNy+1
173			DO i=0,sNx+1
174			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
175			& *recip_hFacC(i,j,k,bi,bj)
176			ENDDO
177	jmc	1.3	ENDDO
178			ELSE
179	jmc	1.9	DO j=0,sNy+1
180			DO i=0,sNx+1
181			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
182			& recip_hFacC(i,j,k,bi,bj)recip_rA(i,j,bi,bj)
183			& Shap_TrLengthShap_TrLength
184			ENDDO
185	jmc	1.3	ENDDO
186			ENDIF
187	jmc	1.9	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
188	jmc	1.3
189	jmc	1.9	C end k,bi,bj loop:
190	jmc	1.3	ENDDO
191			ENDDO
192			ENDDO
193	jmc	1.9	C end loop n=1,nShapTr
194	adcroft	1.2	ENDDO
195
196	jmc	1.4	C F <- [1 - (d_xx+d_yy)^n *deltaT/tau].F
197	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
198			DO bi=myBxLo(myThid),myBxHi(myThid)
199	jmc	1.9	DO k=1,kSize
200			DO j=1,sNy
201			DO i=1,sNx
202	jmc	1.3	field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
203	jmc	1.7	& -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
204			tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau
205	adcroft	1.2	ENDDO
206			ENDDO
207			ENDDO
208			ENDDO
209			ENDDO
210
211	jmc	1.9	c IF (kSize.EQ.Nr) THEN
212			c _EXCH_XYZ_R8( field, myThid )
213			c ELSEIF (kSize.EQ.1) THEN
214			c _EXCH_XY_R8( field, myThid )
215			c ELSE
216			c STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
217			c ENDIF
218	adcroft	1.2
219			ENDIF
220			#endif /* ALLOW_SHAP_FILT */
221
222			RETURN
223			END