/[MITgcm]/MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F
ViewVC logotype

Diff of /MITgcm/pkg/shap_filt/shap_filt_tracer_s2.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by adcroft, Tue May 29 14:01:40 2001 UTC revision 1.11 by jmc, Sun Mar 26 22:56:43 2006 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "SHAP_FILT_OPTIONS.h"  #include "SHAP_FILT_OPTIONS.h"
5    
6        SUBROUTINE SHAP_FILT_TRACER_S2(  CBOP
7       U           field,  C     !ROUTINE: SHAP_FILT_TRACER_S2
8       I           myTime, myThid )  C     !INTERFACE:
9  C     /==========================================================\        SUBROUTINE SHAP_FILT_TRACER_S2(
10  C     | S/R SHAP_FILT_TRACER_S2                                  |       U           field, tmpFld,
11  C     | Applies Shapiro filter to tracer field over one XY slice |       I           nShapTr, exchInOut, kSize,
12  C     | of one tile at a time.                                   |       I           myTime, myIter, myThid )
13  C     \==========================================================/  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    
23    C     !USES:
24        IMPLICIT NONE        IMPLICIT NONE
25    
26  C     == Global variables ===  C     == Global variables ===
# Line 19  C     == Global variables === Line 29  C     == Global variables ===
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "GRID.h"  #include "GRID.h"
31  #include "SHAP_FILT.h"  #include "SHAP_FILT.h"
 #include "SHAP_FILT_TRACER.h"  
32    
33    C     !INPUT/OUTPUT PARAMETERS:
34  C     == Routine arguments  C     == Routine arguments
35        _RL field(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  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          _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          INTEGER nShapTr, exchInOut, kSize
51        _RL     myTime        _RL     myTime
52          INTEGER myIter
53        INTEGER myThid        INTEGER myThid
54    
55  #ifdef ALLOW_SHAP_FILT  #ifdef ALLOW_SHAP_FILT
56    
57    C     !LOCAL VARIABLES:
58  C     == Local variables ==  C     == Local variables ==
59        INTEGER bi,bj,K,I,J,N        INTEGER nShapComput
60          INTEGER bi,bj,k,i,j,n
61        _RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62          _RL tmpFdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63          _RL tmpFdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64    CEOP
65    
66          IF ( exchInOut.LT.0 .OR. exchInOut.GT.3 ) THEN
67             STOP 'S/R SHAP_FILT_TRACER_S2: exchInOut is wrong'
68          ENDIF
69    
70        IF (nShapT.gt.0) THEN        IF (nShapTr.GT.0) THEN
71    C-------
72    C  Apply computational filter ^(nShap-nShapPhys) without grid factor
73    C  then apply Physical filter ^nShapPhys  with grid factors
74    C-------
75            nShapComput = nShapTr - nShapTrPhys
76    
77          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
78           DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
79            DO K=1,Nr            DO k=1,kSize
80             DO J=1,sNy             DO j=1-Oly,sNy+Oly
81              DO I=1,sNx              DO i=1-Olx,sNx+Olx
82               tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)               tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
83              ENDDO              ENDDO
84             ENDDO             ENDDO
# Line 47  C     == Local variables == Line 87  C     == Local variables ==
87          ENDDO          ENDDO
88    
89    
90  C      ( d_xx +d_yy )^n tmpFld  C      ( d_xx +d_yy )^n tmpFld
91    
92         DO N=1,nShapT         DO n=1,nShapTr
93    
94          _EXCH_XYZ_R8( tmpFld, myThid )          IF ( ( MOD(n,2).EQ.1 .OR. Shap_alwaysExchTr ) .AND.
95         &       ( n.GE.2 .OR. MOD(exchInOut,2).EQ.1 )  ) THEN
96             IF (kSize.EQ.Nr) THEN
97              _EXCH_XYZ_R8( tmpFld, myThid )
98             ELSEIF (kSize.EQ.1) THEN
99              _EXCH_XY_R8( tmpFld, myThid )
100             ELSE
101              STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
102             ENDIF
103            ENDIF
104    
105          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
106           DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
107            DO K=1,Nr            DO k=1,kSize
108    
109             DO J=1,sNy  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
             DO I=1,sNx  
              tmpGrd(i,j) = -0.125*(  
      &        ( 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  
110    
111             DO J=1,sNy  C--        Calculate gradient in X direction:
112              DO I=1,sNx  #ifndef ALLOW_AUTODIFF_TAMC
113               tmpFld(i,j,k,bi,bj) = tmpGrd(i,j)             IF ( .NOT.Shap_alwaysExchTr
114              ENDDO       &          .AND. useCubedSphereExchange ) THEN
115    C          to compute d/dx(tmpFld), fill corners with appropriate values:
116                 CALL FILL_CS_CORNER_TR_RL( .TRUE.,
117         &                                 tmpFld(1-Olx,1-Oly,k,bi,bj),
118         &                                 bi,bj, myThid )
119               ENDIF
120    #endif
121               IF ( n.LE.nShapComput ) THEN
122    C-         Computational space: del_i
123                 DO j=0,sNy+1
124                  DO i=0,sNx+2
125                   tmpFdx(i,j) =
126         &            ( 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                   tmpFdx(i,j) =
135         &            ( 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    
142    C--        Calculate gradient in Y direction:
143    #ifndef ALLOW_AUTODIFF_TAMC
144               IF ( .NOT.Shap_alwaysExchTr
145         &          .AND. useCubedSphereExchange ) THEN
146    C          to compute d/dy(tmpFld), fill corners with appropriate values:
147                 CALL FILL_CS_CORNER_TR_RL(.FALSE.,
148         &                                 tmpFld(1-Olx,1-Oly,k,bi,bj),
149         &                                 bi,bj, myThid )
150               ENDIF
151    #endif
152               IF ( n.LE.nShapComput ) THEN
153    C-         Computational space: del_j
154                 DO j=0,sNy+2
155                  DO i=0,sNx+1
156                   tmpFdy(i,j) =
157         &            ( 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                   tmpFdy(i,j) =
166         &            ( 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    
173    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             ENDDO             ENDDO
180    
181    C--        Computational space Filter
182               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    C--        Physical space Filter
189               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                 ENDDO
196               ELSE
197                 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_TrLength*Shap_TrLength
202                  ENDDO
203                 ENDDO
204               ENDIF
205    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
206    
207    C      end k,bi,bj loop:
208            ENDDO            ENDDO
209           ENDDO           ENDDO
210          ENDDO          ENDDO
211    C      end loop n=1,nShapTr
212         ENDDO         ENDDO
213    
214  C      F <-  [1-(d_xx+d_yy)^n]F  C      F <-  [1 - (d_xx+d_yy)^n *deltaT/tau].F
215         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
216          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
217           DO K=1,Nr           DO k=1,kSize
218            DO J=1,sNy            DO j=1,sNy
219             DO I=1,sNx             DO i=1,sNx
220              field(i,j,k,bi,bj)=field(i,j,k,bi,bj)-tmpFld(i,j,k,bi,bj)              field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
221         &            -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             ENDDO             ENDDO
224            ENDDO            ENDDO
225           ENDDO           ENDDO
226          ENDDO          ENDDO
227         ENDDO         ENDDO
228    
229         _EXCH_XYZ_R8( field, myThid )         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    
239        ENDIF        ENDIF
240  #endif /* ALLOW_SHAP_FILT */  #endif /* ALLOW_SHAP_FILT */
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22