pkg/shap_filt/shap_filt_tracer_s2.F

C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.9 2005/10/07 00:27:53 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:
#ifndef ALLOW_AUTODIFF_TAMC
           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
#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:
#ifndef ALLOW_AUTODIFF_TAMC
           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
#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.10	C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.9 2005/10/07 00:27:53 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	jmc	1.10	#ifndef ALLOW_AUTODIFF_TAMC
99	jmc	1.9	IF ( .NOT.Shap_alwaysExchTr
100			& .AND. useCubedSphereExchange ) THEN
101			C to compute d/dx(tmpFld), fill corners with appropriate values:
102			CALL FILL_CS_CORNER_TR_RL( .TRUE.,
103			& tmpFld(1-Olx,1-Oly,k,bi,bj),
104			& bi,bj, myThid )
105			ENDIF
106	jmc	1.10	#endif
107	jmc	1.9	IF ( n.LE.nShapComput ) THEN
108			C- Computational space: del_i
109			DO j=0,sNy+1
110			DO i=0,sNx+2
111			tmpFdx(i,j) =
112			& ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
113			& *_maskW(i,j,k,bi,bj)
114			ENDDO
115			ENDDO
116			ELSE
117			C- Physical space: grad_x
118			DO j=0,sNy+1
119			DO i=0,sNx+2
120			tmpFdx(i,j) =
121			& ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
122			& *_hFacW(i,j,k,bi,bj)
123			& dyG(i,j,bi,bj)recip_dxC(i,j,bi,bj)
124			ENDDO
125			ENDDO
126			ENDIF
127	adcroft	1.2
128	jmc	1.9	C-- Calculate gradient in Y direction:
129	jmc	1.10	#ifndef ALLOW_AUTODIFF_TAMC
130	jmc	1.9	IF ( .NOT.Shap_alwaysExchTr
131			& .AND. useCubedSphereExchange ) THEN
132			C to compute d/dy(tmpFld), fill corners with appropriate values:
133			CALL FILL_CS_CORNER_TR_RL(.FALSE.,
134			& tmpFld(1-Olx,1-Oly,k,bi,bj),
135			& bi,bj, myThid )
136			ENDIF
137	jmc	1.10	#endif
138	jmc	1.9	IF ( n.LE.nShapComput ) THEN
139			C- Computational space: del_j
140			DO j=0,sNy+2
141			DO i=0,sNx+1
142			tmpFdy(i,j) =
143			& ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
144			& *_maskS(i,j,k,bi,bj)
145			ENDDO
146			ENDDO
147			ELSE
148			C- Physical space: grad_y
149			DO j=0,sNy+2
150			DO i=0,sNx+1
151			tmpFdy(i,j) =
152			& ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
153			& *_hFacS(i,j,k,bi,bj)
154			& dxG(i,j,bi,bj)recip_dyC(i,j,bi,bj)
155			ENDDO
156			ENDDO
157			ENDIF
158	jmc	1.3
159	jmc	1.9	C-- Calculate (d_xx + d_yy) tmpFld :
160			DO j=0,sNy+1
161			DO i=0,sNx+1
162			tmpGrd(i,j) = ( tmpFdx(i+1,j) - tmpFdx(i,j) )
163			& + ( tmpFdy(i,j+1) - tmpFdy(i,j) )
164			ENDDO
165	jmc	1.3	ENDDO
166
167	jmc	1.9	C-- Computational space Filter
168			IF ( n.LE.nShapComput ) THEN
169			DO j=0,sNy+1
170			DO i=0,sNx+1
171			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
172			ENDDO
173			ENDDO
174			C-- Physical space Filter
175			ELSEIF (Shap_TrLength.LE.0.) THEN
176			DO j=0,sNy+1
177			DO i=0,sNx+1
178			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
179			& *recip_hFacC(i,j,k,bi,bj)
180			ENDDO
181	jmc	1.3	ENDDO
182			ELSE
183	jmc	1.9	DO j=0,sNy+1
184			DO i=0,sNx+1
185			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
186			& recip_hFacC(i,j,k,bi,bj)recip_rA(i,j,bi,bj)
187			& Shap_TrLengthShap_TrLength
188			ENDDO
189	jmc	1.3	ENDDO
190			ENDIF
191	jmc	1.9	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
192	jmc	1.3
193	jmc	1.9	C end k,bi,bj loop:
194	jmc	1.3	ENDDO
195			ENDDO
196			ENDDO
197	jmc	1.9	C end loop n=1,nShapTr
198	adcroft	1.2	ENDDO
199
200	jmc	1.4	C F <- [1 - (d_xx+d_yy)^n *deltaT/tau].F
201	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
202			DO bi=myBxLo(myThid),myBxHi(myThid)
203	jmc	1.9	DO k=1,kSize
204			DO j=1,sNy
205			DO i=1,sNx
206	jmc	1.3	field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
207	jmc	1.7	& -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
208			tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau
209	adcroft	1.2	ENDDO
210			ENDDO
211			ENDDO
212			ENDDO
213			ENDDO
214
215	jmc	1.9	c IF (kSize.EQ.Nr) THEN
216			c _EXCH_XYZ_R8( field, myThid )
217			c ELSEIF (kSize.EQ.1) THEN
218			c _EXCH_XY_R8( field, myThid )
219			c ELSE
220			c STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
221			c ENDIF
222	adcroft	1.2
223			ENDIF
224			#endif /* ALLOW_SHAP_FILT */
225
226			RETURN
227			END