/[MITgcm]/MITgcm/pkg/gmredi/gmredi_slope_limit.F
ViewVC logotype

Diff of /MITgcm/pkg/gmredi/gmredi_slope_limit.F

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

revision 1.12 by heimbach, Sun Mar 24 02:33:16 2002 UTC revision 1.13 by heimbach, Thu Nov 14 22:43:49 2002 UTC
# Line 4  Line 4 
4  CStartOfInterface  CStartOfInterface
5        SUBROUTINE GMREDI_SLOPE_LIMIT(        SUBROUTINE GMREDI_SLOPE_LIMIT(
6       I             dSigmaDrReal,       I             dSigmaDrReal,
7       I             depthZ,       I             depthZ,K,
8       U             SlopeX, SlopeY,       U             SlopeX, SlopeY,
9         U             dSigmaDx, dSigmaDy,
10       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
11       I             bi,bj, myThid )       I             bi,bj, myThid )
12  C     /==========================================================\  C     /==========================================================\
# Line 40  C     == Routine arguments == Line 41  C     == Routine arguments ==
41  C  C
42        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44          _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45          _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49        _RL depthZ        _RL depthZ
50        INTEGER bi,bj,myThid        INTEGER bi,bj,K,myThid
51  CEndOfInterface  CEndOfInterface
52    
53  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
# Line 86  C     == Local variables == Line 89  C     == Local variables ==
89        IF (GM_taper_scheme.EQ.'orig' .OR.        IF (GM_taper_scheme.EQ.'orig' .OR.
90       &    GM_taper_scheme.EQ.'clipping') THEN       &    GM_taper_scheme.EQ.'clipping') THEN
91    
92    #if (defined (GM_TAPER_ORIG_CLIPPING))
93    
94  C-      Original implementation in mitgcmuv  C-      Original implementation in mitgcmuv
95  C       (this turns out to be the same as Cox slope clipping)  C       (this turns out to be the same as Cox slope clipping)
96    
97  C-      Cox 1987 "Slope clipping"  C-      Cox 1987 "Slope clipping"
98          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
99           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
100            IF ( SlopeX(i,j)*SlopeX(i,j) + SlopeY(i,j)*SlopeY(i,j)            IF ( DSigmaDX(i,j)*DSigmaDX(i,j) +
101         &        DSigmaDY(i,j)*DSigmaDY(i,j)
102       &        .EQ. 0. ) THEN       &        .EQ. 0. ) THEN
103             gradSmod(i,j) = 0. _d 0             gradSmod(i,j) = 0. _d 0
104            ELSE            ELSE
105             gradSmod(i,j) = sqrt(SlopeX(i,j)*SlopeX(i,j)             gradSmod(i,j) = sqrt(DSigmaDX(i,j)*DSigmaDX(i,j)
106       &                         +SlopeY(i,j)*SlopeY(i,j))       &                         +DSigmaDY(i,j)*DSigmaDY(i,j))
107            ENDIF            ENDIF
108           ENDDO           ENDDO
109          ENDDO          ENDDO
110    
111  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
112  CADJ STORE gradSmod(:,:)     = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE gradSmod(:,:)     = comlev1_bibj, key=ikey, byte=isbyte
113  CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
114  #endif  #endif
115    
116          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
117           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
118            IF (gradSmod(i,j) .NE. 0.) THEN            IF (gradSmod(i,j) .NE. 0.) THEN
119             dSigmaDrLtd(i,j) = -gradSmod(i,j)*GM_rMaxSlope             dSigmaDrLtd(i,j) = -gradSmod(i,j)*GM_rMaxSlope
120             IF (dSigmaDrReal(i,j) .GE. dSigmaDrLtd(i,j))             IF ( dSigmaDrReal(i,j) .GE.
121       &         dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)       &          GM_adjointRescale*dSigmaDrLtd(i,j) )
122         &          dSigmaDrReal(i,j) =
123         &          dSigmaDrLtd(i,j)*GM_adjointRescale
124    ctest           dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
125            ENDIF            ENDIF
126           ENDDO           ENDDO
127          ENDDO          ENDDO
128    
129  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
130  CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
131  CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
132  CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
133  #endif  #endif
134    
135          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
# Line 129  CADJ STORE dSigmaDrReal(:,:) = comlev1_b Line 138  CADJ STORE dSigmaDrReal(:,:) = comlev1_b
138             SlopeX(i,j) = 0. _d 0             SlopeX(i,j) = 0. _d 0
139             SlopeY(i,j) = 0. _d 0             SlopeY(i,j) = 0. _d 0
140            ELSE            ELSE
141             dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)             dRdSigmaLtd(i,j) = 1./( dSigmaDrReal(i,j) )
142             SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)             SlopeX(i,j)=-DSigmaDX(i,j)*dRdSigmaLtd(i,j)
143             SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)             SlopeY(i,j)=-DSigmaDY(i,j)*dRdSigmaLtd(i,j)
144            ENDIF            ENDIF
145           ENDDO           ENDDO
146          ENDDO          ENDDO
147    
148  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
149  CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
150  CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
151  #endif  #endif
152    
153          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
# Line 149  CADJ STORE slopeY(:,:)       = comlev1_b Line 158  CADJ STORE slopeY(:,:)       = comlev1_b
158           ENDDO           ENDDO
159          ENDDO          ENDDO
160    
161    #else
162    
163            STOP 'Need to set runtime flag GM_taper_scheme correctly'
164    
165    #endif /* GM_TAPER_ORIG_CLIPPING */
166    
167          ELSE IF (GM_taper_scheme.EQ.'ac02') THEN
168    
169    #if (defined (GM_TAPER_AC02))
170    
171            DO j=1-Oly+1,sNy+Oly-1
172             DO i=1-Olx+1,sNx+Olx-1
173              dRdSigmaLtd(i,j)=
174         &                      dSigmaDx(i,j)*dSigmaDx(i,j)
175         &                      + dSigmaDy(i,j)*dSigmaDy(i,j)
176         &                      + dSigmaDrReal(i,j)**2
177    ctest
178              IF (dRdSigmaLtd(i,j).NE.0.) then
179                 dRdSigmaLtd(i,j)=1. _d 0
180         &            / ( dRdSigmaLtd(i,j) )
181                 SlopeSqr(i,j)=(dSigmaDx(i,j)*dSigmaDx(i,j)
182         &            +dSigmaDy(i,j)*dSigmaDy(i,j))*dRdSigmaLtd(i,j)
183                 SlopeX(i,j)=-dSigmaDx(i,j)
184         &            *dRdSigmaLtd(i,j)*dSigmaDrReal(i,j)
185                 SlopeY(i,j)=-dSigmaDy(i,j)
186         &            *dRdSigmaLtd(i,j)*dSigmaDrReal(i,j)
187    cph             T11(i,j)=(dSigmaDrReal(i,j)**2)*dRdSigmaLtd(i,j)
188                 taperFct(i,j) = 1. _d 0
189    cph          ELSE
190    cph             SlopeSqr(i,j) = 0. _d 0
191    cph             SlopeX(i,j)   = 0. _d 0
192    cph             SlopeY(i,j)   = 0. _d 0
193    cph             T11(i,j)      = 0. _d 0
194    cph             taperFct(i,j) = 0. _d 0
195              ENDIF
196    cph          IF ( SlopeSqr(i,j) .GT. GM_maxSlope**2 ) THEN
197    cph           taperFct(i,j) = GM_maxSlope**2/SlopeSqr(i,j)
198    cph          ENDIF
199             ENDDO
200            ENDDO
201    
202    #else
203          
204            STOP 'Need to set runtime flag GM_taper_scheme correctly'
205    
206    #endif /* GM_TAPER_AC02 */
207    
208        ELSE        ELSE
209    
210    #if (defined (GM_TAPER_REST))
211    
212  C----------------------------------------------------------------------  C----------------------------------------------------------------------
213    
214  C- Compute the slope, no clipping, but avoid reverse slope in negatively  C- Compute the slope, no clipping, but avoid reverse slope in negatively
215  C                                  stratified (Sigma_Z > 0) region :  C                                  stratified (Sigma_Z > 0) region :
216    
217  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
218  CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
219  #endif  #endif
220    
221          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
# Line 170  CADJ STORE dSigmaDrReal(:,:) = comlev1_b Line 228  CADJ STORE dSigmaDrReal(:,:) = comlev1_b
228          ENDDO          ENDDO
229    
230  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
231  CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE dSigmaDx(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
232  CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE dSigmaDy(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
233  CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
234  #endif  #endif
235    
236          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
237           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
238            IF ( dSigmaDrReal(i,j) .EQ. 0. ) THEN            IF ( dSigmaDrReal(i,j) .EQ. 0. ) THEN
239             IF ( SlopeX(i,j) .NE. 0. )             IF ( dSigmaDx(i,j) .NE. 0. )
240       &      SlopeX(i,j) = SIGN(Small_taper,SlopeX(i,j))       &      SlopeX(i,j) = SIGN(Small_taper,dSigmaDx(i,j))
241             IF ( SlopeY(i,j) .NE. 0. )             IF ( dSigmaDy(i,j) .NE. 0. )
242       &      SlopeY(i,j) = SIGN(Small_taper,SlopeY(i,j))       &      SlopeY(i,j) = SIGN(Small_taper,dSigmaDy(i,j))
243            ELSE            ELSE
244             dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)             dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)
245             SlopeX(i,j) = -SlopeX(i,j)*dRdSigmaLtd(i,j)             SlopeX(i,j) = -dSigmaDx(i,j)*dRdSigmaLtd(i,j)
246             SlopeY(i,j) = -SlopeY(i,j)*dRdSigmaLtd(i,j)             SlopeY(i,j) = -dSigmaDy(i,j)*dRdSigmaLtd(i,j)
247            ENDIF            ENDIF
248           ENDDO           ENDDO
249          ENDDO          ENDDO
250    
251  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
252  CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
253  CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte  cnostore CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
254  #endif  #endif
255    
256          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
# Line 280  C-      Large, Danabasoglu and Doney, JP Line 338  C-      Large, Danabasoglu and Doney, JP
338          STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'          STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'
339         ENDIF         ENDIF
340    
341    #endif
342    
343        ENDIF        ENDIF
344    
345    
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22