pkg/shap_filt/shap_filt_tracer_s2.F

C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.14 2008/10/22 00:28:47 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, exchInOut, kSize,
     I           myTime, myIter, 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     exchInOut :: apply Exchanges to fill overlap region:
C            = 0 : do not apply Exch on neither input nor output field
C            = 1 : apply Exch on input field
C                  (needed if input field has invalid overlap)
C            = 2 : apply Exch on output field (after the filter)
C            = 3 : apply Exch on both input & output field
C     kSize     :: length of 3rd Dim : either =1 (2D field) or =Nr (3D field)
C     myTime    :: Current time in simulation
C     myIter    :: Current iteration number in simulation
C     myThid    :: Thread number for this instance of SHAP_FILT_TRACER_S2
      INTEGER 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)
      INTEGER nShapTr, exchInOut
      _RL     myTime
      INTEGER myIter
      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 ( exchInOut.LT.0 .OR. exchInOut.GT.3 ) THEN
         STOP 'S/R SHAP_FILT_TRACER_S2: exchInOut is wrong'
      ENDIF

      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 ) .AND.
     &       ( n.GE.2 .OR. MOD(exchInOut,2).EQ.1 )  ) THEN
         IF (kSize.EQ.Nr) THEN
          _EXCH_XYZ_RL( tmpFld, myThid )
         ELSEIF (kSize.EQ.1) THEN
          _EXCH_XY_RL( 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( 1, .FALSE.,
     &                                 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( 2, .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

       IF ( exchInOut.GE.2 ) THEN
        IF (kSize.EQ.Nr) THEN
          _EXCH_XYZ_RL( field, myThid )
        ELSEIF (kSize.EQ.1) THEN
          _EXCH_XY_RL( field, myThid )
        ELSE
          STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
        ENDIF
       ENDIF

      ENDIF
#endif /* ALLOW_SHAP_FILT */

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