pkg/shap_filt/shap_filt_tracer_s2.F

C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.7 2005/02/15 00:21:12 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, nShapPhysic
      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 = MAX( 0, nShapTr - nShapTrPhys )
        nShapPhysic = nShapTr - nShapComput

        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,nShapPhysic
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.LE.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
            tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/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.8	C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.7 2005/02/15 00:21:12 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.8	INTEGER nShapComput, nShapPhysic
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.8	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.8	nShapComput = MAX( 0, nShapTr - nShapTrPhys )
61			nShapPhysic = nShapTr - nShapComput
62	adcroft	1.2
63			DO bj=myByLo(myThid),myByHi(myThid)
64			DO bi=myBxLo(myThid),myBxHi(myThid)
65	jmc	1.4	DO K=1,kSize
66	jmc	1.3	DO J=1-Oly,sNy+Oly
67			DO I=1-Olx,sNx+Olx
68	adcroft	1.2	tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
69			ENDDO
70			ENDDO
71			ENDDO
72			ENDDO
73			ENDDO
74
75
76			C ( d_xx +d_yy )^n tmpFld
77
78	jmc	1.3	C-- Computational Filter
79			DO N=1,nShapComput
80	adcroft	1.2
81	jmc	1.4	IF (kSize.EQ.Nr) THEN
82			_EXCH_XYZ_R8( tmpFld, myThid )
83			ELSE
84			_EXCH_XY_R8( tmpFld, myThid )
85			ENDIF
86	jmc	1.8
87	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
88			DO bi=myBxLo(myThid),myBxHi(myThid)
89	jmc	1.4	DO K=1,kSize
90	jmc	1.8
91	adcroft	1.2	DO J=1,sNy
92			DO I=1,sNx
93	jmc	1.3	tmpGrd(i,j) =
94	adcroft	1.2	& ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )
95			& *_maskW(i+1,j,k,bi,bj)
96			& -( tmpFld( i ,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
97	jmc	1.3	& *_maskW( i ,j,k,bi,bj)
98	adcroft	1.2	& +( tmpFld(i,j+1,k,bi,bj)-tmpFld(i, j ,k,bi,bj) )
99			& *_maskS(i,j+1,k,bi,bj)
100			& -( tmpFld(i, j ,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
101	jmc	1.3	& *_maskS(i, j ,k,bi,bj)
102	adcroft	1.2	ENDDO
103			ENDDO
104
105			DO J=1,sNy
106			DO I=1,sNx
107	jmc	1.3	tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
108	adcroft	1.2	ENDDO
109			ENDDO
110
111			ENDDO
112			ENDDO
113			ENDDO
114	jmc	1.3	C end loop N=1,nShapComput
115			ENDDO
116
117	jmc	1.8	DO N=1,nShapPhysic
118	jmc	1.3	C-- Physical space Filter
119
120	jmc	1.4	IF (kSize.EQ.Nr) THEN
121			_EXCH_XYZ_R8( tmpFld, myThid )
122			ELSE
123			_EXCH_XY_R8( tmpFld, myThid )
124			ENDIF
125	adcroft	1.2
126	jmc	1.3	DO bj=myByLo(myThid),myByHi(myThid)
127			DO bi=myBxLo(myThid),myBxHi(myThid)
128	jmc	1.4	DO K=1,kSize
129	jmc	1.3
130			DO J=1,sNy
131			DO I=1,sNx
132			tmpGrd(i,j) =
133			& ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )
134			& *_hFacW(i+1,j,k,bi,bj)
135			& *DYG(i+1,j,bi,bj)
136			& *recip_DXC(i+1,j,bi,bj)
137			& -( tmpFld( i ,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
138			& *_hFacW( i ,j,k,bi,bj)
139			& *DYG( i ,j,bi,bj)
140			& *recip_DXC( i ,j,bi,bj)
141			& +( tmpFld(i,j+1,k,bi,bj)-tmpFld(i, j ,k,bi,bj) )
142			& *_hFacS(i,j+1,k,bi,bj)
143			& *DXG(i,j+1,bi,bj)
144			& *recip_DYC(i,j+1,bi,bj)
145			& -( tmpFld(i, j ,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
146			& *_hFacS(i, j ,k,bi,bj)
147			& *DXG(i, j ,bi,bj)
148			& *recip_DYC(i, j ,bi,bj)
149			ENDDO
150			ENDDO
151
152	jmc	1.8	IF (Shap_TrLength.LE.0.) THEN
153	jmc	1.3	DO J=1,sNy
154			DO I=1,sNx
155			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
156			& *recip_hFacC(i,j,k,bi,bj)
157			ENDDO
158			ENDDO
159			ELSE
160			DO J=1,sNy
161			DO I=1,sNx
162			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
163			& recip_hFacC(i,j,k,bi,bj)recip_rA(i,j,bi,bj)
164			& Shap_TrLengthShap_TrLength
165			ENDDO
166			ENDDO
167			ENDIF
168
169			ENDDO
170			ENDDO
171			ENDDO
172			C end loop N=1,nShapTrPhys
173	adcroft	1.2	ENDDO
174
175	jmc	1.4	C F <- [1 - (d_xx+d_yy)^n *deltaT/tau].F
176	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
177			DO bi=myBxLo(myThid),myBxHi(myThid)
178	jmc	1.4	DO K=1,kSize
179	adcroft	1.2	DO J=1,sNy
180			DO I=1,sNx
181	jmc	1.3	field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
182	jmc	1.7	& -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
183			tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau
184	adcroft	1.2	ENDDO
185			ENDDO
186			ENDDO
187			ENDDO
188			ENDDO
189
190	jmc	1.4	IF (kSize.EQ.Nr) THEN
191			_EXCH_XYZ_R8( field, myThid )
192			ELSEIF (kSize.EQ.1) THEN
193			_EXCH_XY_R8( field, myThid )
194			ELSE
195			STOP 'S/R SHAP_FILT_TRACER_S4: kSize is wrong'
196			ENDIF
197	adcroft	1.2
198			ENDIF
199			#endif /* ALLOW_SHAP_FILT */
200
201			RETURN
202			END