/[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.8 by jmc, Tue Sep 27 22:11:06 2005 UTC revision 1.9 by jmc, Fri Oct 7 00:27:53 2005 UTC
# Line 47  C     myThid :: Thread number for this i Line 47  C     myThid :: Thread number for this i
47    
48  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
49  C     == Local variables ==  C     == Local variables ==
50        INTEGER nShapComput, nShapPhysic        INTEGER nShapComput
51        INTEGER bi,bj,K,I,J,N        INTEGER bi,bj,k,i,j,n
52        _RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tmpGrd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53          _RL tmpFdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54          _RL tmpFdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55  CEOP  CEOP
56    
57        IF (nShapTr.GT.0) THEN        IF (nShapTr.GT.0) THEN
58  C-------  C-------
59  C  Apply computational filter ^(nShap-nShapPhys) without grid factor  C  Apply computational filter ^(nShap-nShapPhys) without grid factor
60  C  then apply Physical filter ^nShapPhys  with grid factors  C  then apply Physical filter ^nShapPhys  with grid factors
61  C-------  C-------
62          nShapComput = MAX( 0, nShapTr - nShapTrPhys )          nShapComput = nShapTr - nShapTrPhys
         nShapPhysic = nShapTr - nShapComput  
63    
64          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
65           DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
66            DO K=1,kSize            DO k=1,kSize
67             DO J=1-Oly,sNy+Oly             DO j=1-Oly,sNy+Oly
68              DO I=1-Olx,sNx+Olx              DO i=1-Olx,sNx+Olx
69               tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)               tmpFld(i,j,k,bi,bj)=field(i,j,k,bi,bj)
70              ENDDO              ENDDO
71             ENDDO             ENDDO
# Line 73  C------- Line 74  C-------
74          ENDDO          ENDDO
75    
76    
77  C      ( d_xx +d_yy )^n tmpFld  C      ( d_xx +d_yy )^n tmpFld
78    
79  C-- Computational Filter         DO n=1,nShapTr
        DO N=1,nShapComput  
80    
81          IF (kSize.EQ.Nr) THEN          IF ( MOD(n,2).EQ.1 .OR. Shap_alwaysExchTr ) THEN
82             IF (kSize.EQ.Nr) THEN
83            _EXCH_XYZ_R8( tmpFld, myThid )            _EXCH_XYZ_R8( tmpFld, myThid )
84          ELSE           ELSEIF (kSize.EQ.1) THEN
85            _EXCH_XY_R8( tmpFld, myThid )            _EXCH_XY_R8( tmpFld, myThid )
86             ELSE
87              STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
88             ENDIF
89          ENDIF          ENDIF
90    
91          DO bj=myByLo(myThid),myByHi(myThid)          DO bj=myByLo(myThid),myByHi(myThid)
92           DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
93            DO K=1,kSize            DO k=1,kSize
94    
95             DO J=1,sNy  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
             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  
96    
97         DO N=1,nShapPhysic  C--        Calculate gradient in X direction:
98  C-- Physical space Filter             IF ( .NOT.Shap_alwaysExchTr
99         &          .AND. useCubedSphereExchange ) THEN
100          IF (kSize.EQ.Nr) THEN  C          to compute d/dx(tmpFld), fill corners with appropriate values:
101            _EXCH_XYZ_R8( tmpFld, myThid )               CALL FILL_CS_CORNER_TR_RL( .TRUE.,
102          ELSE       &                                 tmpFld(1-Olx,1-Oly,k,bi,bj),
103            _EXCH_XY_R8( tmpFld, myThid )       &                                 bi,bj, myThid )
104          ENDIF             ENDIF
105               IF ( n.LE.nShapComput ) THEN
106    C-         Computational space: del_i
107                 DO j=0,sNy+1
108                  DO i=0,sNx+2
109                   tmpFdx(i,j) =
110         &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
111         &                     *_maskW(i,j,k,bi,bj)
112                  ENDDO
113                 ENDDO
114               ELSE
115    C-         Physical space: grad_x
116                 DO j=0,sNy+1
117                  DO i=0,sNx+2
118                   tmpFdx(i,j) =
119         &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i-1,j,k,bi,bj) )
120         &                     *_hFacW(i,j,k,bi,bj)
121         &                     *dyG(i,j,bi,bj)*recip_dxC(i,j,bi,bj)
122                  ENDDO
123                 ENDDO
124               ENDIF
125    
126          DO bj=myByLo(myThid),myByHi(myThid)  C--        Calculate gradient in Y direction:
127           DO bi=myBxLo(myThid),myBxHi(myThid)             IF ( .NOT.Shap_alwaysExchTr
128            DO K=1,kSize       &          .AND. useCubedSphereExchange ) THEN
129    C          to compute d/dy(tmpFld), fill corners with appropriate values:
130                 CALL FILL_CS_CORNER_TR_RL(.FALSE.,
131         &                                 tmpFld(1-Olx,1-Oly,k,bi,bj),
132         &                                 bi,bj, myThid )
133               ENDIF
134               IF ( n.LE.nShapComput ) THEN
135    C-         Computational space: del_j
136                 DO j=0,sNy+2
137                  DO i=0,sNx+1
138                   tmpFdy(i,j) =
139         &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
140         &                     *_maskS(i,j,k,bi,bj)
141                  ENDDO
142                 ENDDO
143               ELSE
144    C-         Physical space: grad_y
145                 DO j=0,sNy+2
146                  DO i=0,sNx+1
147                   tmpFdy(i,j) =
148         &            ( tmpFld(i,j,k,bi,bj)-tmpFld(i,j-1,k,bi,bj) )
149         &                     *_hFacS(i,j,k,bi,bj)
150         &                     *dxG(i,j,bi,bj)*recip_dyC(i,j,bi,bj)
151                  ENDDO
152                 ENDDO
153               ENDIF
154    
155             DO J=1,sNy  C--        Calculate (d_xx + d_yy) tmpFld :
156              DO I=1,sNx             DO j=0,sNy+1
157               tmpGrd(i,j) =               DO i=0,sNx+1
158       &        ( tmpFld(i+1,j,k,bi,bj)-tmpFld( i ,j,k,bi,bj) )                 tmpGrd(i,j) = ( tmpFdx(i+1,j) - tmpFdx(i,j) )
159       &            *_hFacW(i+1,j,k,bi,bj)       &                     + ( tmpFdy(i,j+1) - tmpFdy(i,j) )
160       &            *DYG(i+1,j,bi,bj)               ENDDO
      &            *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  
161             ENDDO             ENDDO
162    
163             IF (Shap_TrLength.LE.0.) THEN  C--        Computational space Filter
164              DO J=1,sNy             IF ( n.LE.nShapComput ) THEN
165               DO I=1,sNx               DO j=0,sNy+1
166                tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)                DO i=0,sNx+1
167       &           *recip_hFacC(i,j,k,bi,bj)                 tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
168                  ENDDO
169                 ENDDO
170    C--        Physical space Filter
171               ELSEIF (Shap_TrLength.LE.0.) THEN
172                 DO j=0,sNy+1
173                  DO i=0,sNx+1
174                   tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
175         &             *recip_hFacC(i,j,k,bi,bj)
176                  ENDDO
177               ENDDO               ENDDO
             ENDDO  
178             ELSE             ELSE
179              DO J=1,sNy               DO j=0,sNy+1
180               DO I=1,sNx                DO i=0,sNx+1
181                tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)                 tmpFld(i,j,k,bi,bj) = -0.125*tmpGrd(i,j)
182       &           *recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)       &             *recip_hFacC(i,j,k,bi,bj)*recip_rA(i,j,bi,bj)
183       &           *Shap_TrLength*Shap_TrLength       &             *Shap_TrLength*Shap_TrLength
184                  ENDDO
185               ENDDO               ENDDO
             ENDDO  
186             ENDIF             ENDIF
187    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
188    
189    C      end k,bi,bj loop:
190            ENDDO            ENDDO
191           ENDDO           ENDDO
192          ENDDO          ENDDO
193  C      end loop N=1,nShapTrPhys  C      end loop n=1,nShapTr
194         ENDDO         ENDDO
195    
196  C      F <-  [1 - (d_xx+d_yy)^n *deltaT/tau].F  C      F <-  [1 - (d_xx+d_yy)^n *deltaT/tau].F
197         DO bj=myByLo(myThid),myByHi(myThid)         DO bj=myByLo(myThid),myByHi(myThid)
198          DO bi=myBxLo(myThid),myBxHi(myThid)          DO bi=myBxLo(myThid),myBxHi(myThid)
199           DO K=1,kSize           DO k=1,kSize
200            DO J=1,sNy            DO j=1,sNy
201             DO I=1,sNx             DO i=1,sNx
202              field(i,j,k,bi,bj)=field(i,j,k,bi,bj)              field(i,j,k,bi,bj)=field(i,j,k,bi,bj)
203       &            -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau       &            -tmpFld(i,j,k,bi,bj)*dTtracerLev(1)/Shap_Trtau
204              tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau              tmpFld(i,j,k,bi,bj)= -tmpFld(i,j,k,bi,bj)/Shap_Trtau
# Line 187  C      F <-  [1 - (d_xx+d_yy)^n *deltaT/ Line 208  C      F <-  [1 - (d_xx+d_yy)^n *deltaT/
208          ENDDO          ENDDO
209         ENDDO         ENDDO
210    
211          IF (kSize.EQ.Nr) THEN  c       IF (kSize.EQ.Nr) THEN
212            _EXCH_XYZ_R8( field, myThid )  c         _EXCH_XYZ_R8( field, myThid )
213          ELSEIF (kSize.EQ.1) THEN  c       ELSEIF (kSize.EQ.1) THEN
214            _EXCH_XY_R8( field, myThid )  c         _EXCH_XY_R8( field, myThid )
215          ELSE  c       ELSE
216            STOP 'S/R SHAP_FILT_TRACER_S4: kSize is wrong'  c         STOP 'S/R SHAP_FILT_TRACER_S2: kSize is wrong'
217          ENDIF  c       ENDIF
218    
219        ENDIF        ENDIF
220  #endif /* ALLOW_SHAP_FILT */  #endif /* ALLOW_SHAP_FILT */
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22