pkg/shap_filt/shap_filt_tracer_s2.F

C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.10 2005/10/07 14:19:42 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
      _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, kSize
      _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_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

       IF ( exchInOut.GE.2 ) THEN
        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_S2: kSize is wrong'
        ENDIF
       ENDIF

      ENDIF
#endif /* ALLOW_SHAP_FILT */

      RETURN
      END
1	jmc	1.11	C $Header: /u/gcmpack/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F,v 1.10 2005/10/07 14:19:42 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.4	_RL field(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
49			_RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,kSize,nSx,nSy)
50	jmc	1.11	INTEGER nShapTr, exchInOut, kSize
51	adcroft	1.2	_RL myTime
52	jmc	1.11	INTEGER myIter
53	adcroft	1.2	INTEGER myThid
54	jmc	1.8
55	adcroft	1.2	#ifdef ALLOW_SHAP_FILT
56
57	jmc	1.4	C !LOCAL VARIABLES:
58	adcroft	1.2	C == Local variables ==
59	jmc	1.9	INTEGER nShapComput
60			INTEGER bi,bj,k,i,j,n
61	adcroft	1.2	_RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62	jmc	1.9	_RL tmpFdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63			_RL tmpFdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64	jmc	1.4	CEOP
65	adcroft	1.2
66	jmc	1.11	IF ( exchInOut.LT.0 .OR. exchInOut.GT.3 ) THEN
67			STOP 'S/R SHAP_FILT_TRACER_S2: exchInOut is wrong'
68			ENDIF
69
70	jmc	1.8	IF (nShapTr.GT.0) THEN
71	jmc	1.3	C-------
72			C Apply computational filter ^(nShap-nShapPhys) without grid factor
73	jmc	1.9	C then apply Physical filter ^nShapPhys with grid factors
74	jmc	1.3	C-------
75	jmc	1.9	nShapComput = nShapTr - nShapTrPhys
76	adcroft	1.2
77			DO bj=myByLo(myThid),myByHi(myThid)
78			DO bi=myBxLo(myThid),myBxHi(myThid)
79	jmc	1.9	DO k=1,kSize
80			DO j=1-Oly,sNy+Oly
81			DO i=1-Olx,sNx+Olx
82	adcroft	1.2	tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
83			ENDDO
84			ENDDO
85			ENDDO
86			ENDDO
87			ENDDO
88
89
90	jmc	1.9	C ( d_xx +d_yy )^n tmpFld
91	adcroft	1.2
92	jmc	1.9	DO n=1,nShapTr
93	adcroft	1.2
94	jmc	1.11	IF ( ( MOD(n,2).EQ.1 .OR. Shap_alwaysExchTr ) .AND.
95			& ( n.GE.2 .OR. MOD(exchInOut,2).EQ.1 ) ) THEN
96	jmc	1.9	IF (kSize.EQ.Nr) THEN
97	jmc	1.4	_EXCH_XYZ_R8( tmpFld, myThid )
98	jmc	1.9	ELSEIF (kSize.EQ.1) THEN
99	jmc	1.4	_EXCH_XY_R8( tmpFld, myThid )
100	jmc	1.9	ELSE
101			STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
102			ENDIF
103	jmc	1.4	ENDIF
104	jmc	1.8
105	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
106			DO bi=myBxLo(myThid),myBxHi(myThid)
107	jmc	1.9	DO k=1,kSize
108	jmc	1.8
109	jmc	1.9	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
110	adcroft	1.2
111	jmc	1.9	C-- Calculate gradient in X direction:
112	jmc	1.10	#ifndef ALLOW_AUTODIFF_TAMC
113	jmc	1.9	IF ( .NOT.Shap_alwaysExchTr
114			& .AND. useCubedSphereExchange ) THEN
115			C to compute d/dx(tmpFld), fill corners with appropriate values:
116	jmc	1.11	CALL FILL_CS_CORNER_TR_RL( .TRUE.,
117			& tmpFld(1-Olx,1-Oly,k,bi,bj),
118	jmc	1.9	& bi,bj, myThid )
119			ENDIF
120	jmc	1.10	#endif
121	jmc	1.9	IF ( n.LE.nShapComput ) THEN
122			C- Computational space: del_i
123			DO j=0,sNy+1
124			DO i=0,sNx+2
125	jmc	1.11	tmpFdx(i,j) =
126	jmc	1.9	& ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
127			& *_maskW(i,j,k,bi,bj)
128			ENDDO
129			ENDDO
130			ELSE
131			C- Physical space: grad_x
132			DO j=0,sNy+1
133			DO i=0,sNx+2
134	jmc	1.11	tmpFdx(i,j) =
135	jmc	1.9	& ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
136			& *_hFacW(i,j,k,bi,bj)
137			& dyG(i,j,bi,bj)recip_dxC(i,j,bi,bj)
138			ENDDO
139			ENDDO
140			ENDIF
141	adcroft	1.2
142	jmc	1.9	C-- Calculate gradient in Y direction:
143	jmc	1.10	#ifndef ALLOW_AUTODIFF_TAMC
144	jmc	1.9	IF ( .NOT.Shap_alwaysExchTr
145			& .AND. useCubedSphereExchange ) THEN
146			C to compute d/dy(tmpFld), fill corners with appropriate values:
147	jmc	1.11	CALL FILL_CS_CORNER_TR_RL(.FALSE.,
148			& tmpFld(1-Olx,1-Oly,k,bi,bj),
149	jmc	1.9	& bi,bj, myThid )
150			ENDIF
151	jmc	1.10	#endif
152	jmc	1.9	IF ( n.LE.nShapComput ) THEN
153			C- Computational space: del_j
154			DO j=0,sNy+2
155			DO i=0,sNx+1
156	jmc	1.11	tmpFdy(i,j) =
157	jmc	1.9	& ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
158			& *_maskS(i,j,k,bi,bj)
159			ENDDO
160			ENDDO
161			ELSE
162			C- Physical space: grad_y
163			DO j=0,sNy+2
164			DO i=0,sNx+1
165	jmc	1.11	tmpFdy(i,j) =
166	jmc	1.9	& ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
167			& *_hFacS(i,j,k,bi,bj)
168			& dxG(i,j,bi,bj)recip_dyC(i,j,bi,bj)
169			ENDDO
170			ENDDO
171			ENDIF
172	jmc	1.3
173	jmc	1.9	C-- Calculate (d_xx + d_yy) tmpFld :
174			DO j=0,sNy+1
175			DO i=0,sNx+1
176			tmpGrd(i,j) = ( tmpFdx(i+1,j) - tmpFdx(i,j) )
177			& + ( tmpFdy(i,j+1) - tmpFdy(i,j) )
178			ENDDO
179	jmc	1.3	ENDDO
180
181	jmc	1.11	C-- Computational space Filter
182	jmc	1.9	IF ( n.LE.nShapComput ) THEN
183			DO j=0,sNy+1
184			DO i=0,sNx+1
185			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
186			ENDDO
187			ENDDO
188	jmc	1.11	C-- Physical space Filter
189	jmc	1.9	ELSEIF (Shap_TrLength.LE.0.) THEN
190			DO j=0,sNy+1
191			DO i=0,sNx+1
192			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
193			& *recip_hFacC(i,j,k,bi,bj)
194			ENDDO
195	jmc	1.3	ENDDO
196			ELSE
197	jmc	1.9	DO j=0,sNy+1
198			DO i=0,sNx+1
199			tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
200			& recip_hFacC(i,j,k,bi,bj)recip_rA(i,j,bi,bj)
201			& Shap_TrLengthShap_TrLength
202			ENDDO
203	jmc	1.3	ENDDO
204			ENDIF
205	jmc	1.9	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
206	jmc	1.3
207	jmc	1.9	C end k,bi,bj loop:
208	jmc	1.3	ENDDO
209			ENDDO
210			ENDDO
211	jmc	1.9	C end loop n=1,nShapTr
212	adcroft	1.2	ENDDO
213
214	jmc	1.4	C F <- [1 - (d_xx+d_yy)^n *deltaT/tau].F
215	adcroft	1.2	DO bj=myByLo(myThid),myByHi(myThid)
216			DO bi=myBxLo(myThid),myBxHi(myThid)
217	jmc	1.9	DO k=1,kSize
218			DO j=1,sNy
219			DO i=1,sNx
220	jmc	1.3	field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
221	jmc	1.7	& -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
222			tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau
223	adcroft	1.2	ENDDO
224			ENDDO
225			ENDDO
226			ENDDO
227			ENDDO
228
229	jmc	1.11	IF ( exchInOut.GE.2 ) THEN
230			IF (kSize.EQ.Nr) THEN
231			_EXCH_XYZ_R8( field, myThid )
232			ELSEIF (kSize.EQ.1) THEN
233			_EXCH_XY_R8( field, myThid )
234			ELSE
235			STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
236			ENDIF
237			ENDIF
238	adcroft	1.2
239			ENDIF
240			#endif /* ALLOW_SHAP_FILT */
241
242			RETURN
243			END