pkg/shap_filt/shap_filt_tracer_s2.F

C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.5 2004/06/26 01:15:30 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)
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 

C-- Computational Filter 
       DO N=1,nShapComput

        IF (kSize.EQ.Nr) THEN
          _EXCH_XYZ_R8( tmpFld, myThid )
        ELSE
          _EXCH_XY_R8( tmpFld, myThid )
        ENDIF
 
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO K=1,kSize
 
           DO J=1,sNy
            DO I=1,sNx
             tmpGrd(i,j) =
     &        ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )
     &            *_maskW(i+1,j,k,bi,bj)
     &       -( tmpFld( i ,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
     &            *_maskW( i ,j,k,bi,bj)
     &       +( tmpFld(i,j+1,k,bi,bj)-tmpFld(i, j ,k,bi,bj) )
     &            *_maskS(i,j+1,k,bi,bj)
     &       -( tmpFld(i, j ,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
     &            *_maskS(i, j ,k,bi,bj)
            ENDDO
           ENDDO

           DO J=1,sNy
            DO I=1,sNx
             tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
            ENDDO
           ENDDO

          ENDDO
         ENDDO
        ENDDO
C      end loop N=1,nShapComput
       ENDDO

       DO N=1,nShapTrPhys
C-- Physical space Filter 

        IF (kSize.EQ.Nr) THEN
          _EXCH_XYZ_R8( tmpFld, myThid )
        ELSE
          _EXCH_XY_R8( tmpFld, myThid )
        ENDIF

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

           DO J=1,sNy
            DO I=1,sNx
             tmpGrd(i,j) =
     &        ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )
     &            *_hFacW(i+1,j,k,bi,bj)
     &            *DYG(i+1,j,bi,bj)
     &            *recip_DXC(i+1,j,bi,bj)
     &       -( 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)
     &       +( tmpFld(i,j+1,k,bi,bj)-tmpFld(i, j ,k,bi,bj) )
     &            *_hFacS(i,j+1,k,bi,bj)
     &            *DXG(i,j+1,bi,bj)
     &            *recip_DYC(i,j+1,bi,bj)
     &       -( 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

           IF (Shap_TrLength.EQ.0.) THEN
            DO J=1,sNy
             DO I=1,sNx
              tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
     &           *recip_hFacC(i,j,k,bi,bj)
             ENDDO
            ENDDO
           ELSE
            DO J=1,sNy
             DO I=1,sNx
              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

          ENDDO
         ENDDO
        ENDDO
C      end loop N=1,nShapTrPhys
       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
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO

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

      ENDIF
#endif /* ALLOW_SHAP_FILT */

      RETURN
      END
1	jmc	1.6	C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.5 2004/06/26 01:15:30 jmc Exp $
2	jmc	1.3	C $Name: $
3	adcroft	1.2
4			#include "SHAP_FILT_OPTIONS.h"
5	jmc	1.4
6			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
22			C !USES:
23	adcroft	1.2	IMPLICIT NONE
24	jmc	1.4
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.4
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.3	INTEGER nShapComput
51	adcroft	1.2	INTEGER bi,bj,K,I,J,N
52			_RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53	jmc	1.4	CEOP
54	adcroft	1.2
55	jmc	1.5	IF (nShapTr.gt.0) THEN
56	jmc	1.3	C-------
57			C Apply computational filter ^(nShap-nShapPhys) without grid factor
58			C then apply Physical filter ^nShapPhys with grid factors
59			C-------
60	jmc	1.5	nShapComput = nShapTr - nShapTrPhys
61	adcroft	1.2
62			DO bj=myByLo(myThid),myByHi(myThid)
63			DO bi=myBxLo(myThid),myBxHi(myThid)
64	jmc	1.4	DO K=1,kSize
65	jmc	1.3	DO J=1-Oly,sNy+Oly
66			DO I=1-Olx,sNx+Olx
67	adcroft	1.2	tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
68			ENDDO
69			ENDDO
70			ENDDO
71			ENDDO
72			ENDDO
73
74
75			C ( d_xx +d_yy )^n tmpFld
76
77	jmc	1.3	C-- Computational Filter
78			DO N=1,nShapComput
79	adcroft	1.2
80	jmc	1.4	IF (kSize.EQ.Nr) THEN
81			_EXCH_XYZ_R8( tmpFld, myThid )
82			ELSE
83			_EXCH_XY_R8( tmpFld, myThid )
84			ENDIF
85	jmc	1.3
86	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
87			DO bi=myBxLo(myThid),myBxHi(myThid)
88	jmc	1.4	DO K=1,kSize
89	jmc	1.3
90	adcroft	1.2	DO J=1,sNy
91			DO I=1,sNx
92	jmc	1.3	tmpGrd(i,j) =
93	adcroft	1.2	& ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )
94			& *_maskW(i+1,j,k,bi,bj)
95			& -( tmpFld( i ,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
96	jmc	1.3	& *_maskW( i ,j,k,bi,bj)
97	adcroft	1.2	& +( tmpFld(i,j+1,k,bi,bj)-tmpFld(i, j ,k,bi,bj) )
98			& *_maskS(i,j+1,k,bi,bj)
99			& -( tmpFld(i, j ,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
100	jmc	1.3	& *_maskS(i, j ,k,bi,bj)
101	adcroft	1.2	ENDDO
102			ENDDO
103
104			DO J=1,sNy
105			DO I=1,sNx
106	jmc	1.3	tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
107	adcroft	1.2	ENDDO
108			ENDDO
109
110			ENDDO
111			ENDDO
112			ENDDO
113	jmc	1.3	C end loop N=1,nShapComput
114			ENDDO
115
116			DO N=1,nShapTrPhys
117			C-- Physical space Filter
118
119	jmc	1.4	IF (kSize.EQ.Nr) THEN
120			_EXCH_XYZ_R8( tmpFld, myThid )
121			ELSE
122			_EXCH_XY_R8( tmpFld, myThid )
123			ENDIF
124	adcroft	1.2
125	jmc	1.3	DO bj=myByLo(myThid),myByHi(myThid)
126			DO bi=myBxLo(myThid),myBxHi(myThid)
127	jmc	1.4	DO K=1,kSize
128	jmc	1.3
129			DO J=1,sNy
130			DO I=1,sNx
131			tmpGrd(i,j) =
132			& ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )
133			& *_hFacW(i+1,j,k,bi,bj)
134			& *DYG(i+1,j,bi,bj)
135			& *recip_DXC(i+1,j,bi,bj)
136			& -( tmpFld( i ,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
137			& *_hFacW( i ,j,k,bi,bj)
138			& *DYG( i ,j,bi,bj)
139			& *recip_DXC( i ,j,bi,bj)
140			& +( tmpFld(i,j+1,k,bi,bj)-tmpFld(i, j ,k,bi,bj) )
141			& *_hFacS(i,j+1,k,bi,bj)
142			& *DXG(i,j+1,bi,bj)
143			& *recip_DYC(i,j+1,bi,bj)
144			& -( tmpFld(i, j ,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
145			& *_hFacS(i, j ,k,bi,bj)
146			& *DXG(i, j ,bi,bj)
147			& *recip_DYC(i, j ,bi,bj)
148			ENDDO
149			ENDDO
150
151			IF (Shap_TrLength.EQ.0.) THEN
152			DO J=1,sNy
153			DO I=1,sNx
154			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
155			& *recip_hFacC(i,j,k,bi,bj)
156			ENDDO
157			ENDDO
158			ELSE
159			DO J=1,sNy
160			DO I=1,sNx
161			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
162			& recip_hFacC(i,j,k,bi,bj)recip_rA(i,j,bi,bj)
163			& Shap_TrLengthShap_TrLength
164			ENDDO
165			ENDDO
166			ENDIF
167
168			ENDDO
169			ENDDO
170			ENDDO
171			C end loop N=1,nShapTrPhys
172	adcroft	1.2	ENDDO
173
174	jmc	1.4	C F <- [1 - (d_xx+d_yy)^n *deltaT/tau].F
175	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
176			DO bi=myBxLo(myThid),myBxHi(myThid)
177	jmc	1.4	DO K=1,kSize
178	adcroft	1.2	DO J=1,sNy
179			DO I=1,sNx
180	jmc	1.3	field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
181	jmc	1.6	& -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
182	adcroft	1.2	ENDDO
183			ENDDO
184			ENDDO
185			ENDDO
186			ENDDO
187
188	jmc	1.4	IF (kSize.EQ.Nr) THEN
189			_EXCH_XYZ_R8( field, myThid )
190			ELSEIF (kSize.EQ.1) THEN
191			_EXCH_XY_R8( field, myThid )
192			ELSE
193			STOP 'S/R SHAP_FILT_TRACER_S4: kSize is wrong'
194			ENDIF
195	adcroft	1.2
196			ENDIF
197			#endif /* ALLOW_SHAP_FILT */
198
199			RETURN
200			END