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

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

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

revision 1.1 by jmc, Sun Dec 16 18:54:49 2001 UTC revision 1.12 by jmc, Tue Dec 8 21:42:22 2009 UTC
# Line 4  C $Name$ Line 4  C $Name$
4  #include "GMREDI_OPTIONS.h"  #include "GMREDI_OPTIONS.h"
5    
6  CStartOfInterface  CStartOfInterface
7        SUBROUTINE GMREDI_SLOPE_PSI_B(        SUBROUTINE GMREDI_SLOPE_PSI(
      I             dSigmaDrW,dSigmaDrS,  
      I             depthZ,  
      U             SlopeX, SlopeY,  
8       O             taperX, taperY,       O             taperX, taperY,
9         U             SlopeX, SlopeY,
10         U             dSigmaDrW,dSigmaDrS,
11         I             LrhoW, LrhoS, depthZ, K,
12       I             bi,bj, myThid )       I             bi,bj, myThid )
13  C     /==========================================================\  C     /==========================================================\
14  C     | SUBROUTINE GMREDI_SLOPE_PSI_B                            |  C     | SUBROUTINE GMREDI_SLOPE_PSI                              |
15  C     | o Calculate slopes for use in GM/Redi tensor             |  C     | o Calculate slopes for use in GM/Redi tensor             |
16  C     |==========================================================|  C     |==========================================================|
17  C     | On entry:                                                |  C     | On entry:                                                |
18  C     |            dSigmaDrW conatins the d/dz Sigma          |  C     |            dSigmaDrW,S  contains the d/dz Sigma          |
19  C     |            SlopeX/Y     contains X/Y gradients of sigma  |  C     |            SlopeX/Y     contains X/Y gradients of sigma  |
20  C     |            depthZ       conatins the height (m) of level |  C     |            depthZ       contains the depth (< 0 !) [m]   |
21  C     | On exit:                                                 |  C     | On exit:                                                 |
22  C     |            dSigmaDrW conatins the effective dSig/dz   |  C     |            dSigmaDrW,S  contains the effective dSig/dz   |
23  C     |            SlopeX/Y     contains X/Y slopes              |  C     |            SlopeX/Y     contains X/Y slopes              |
24  C     |            taperFct     contains tapering funct. value ; |  C     |            taperFct     contains tapering funct. value ; |
25  C     |                         = 1 when using no tapering       |  C     |                         = 1 when using no tapering       |
# Line 32  C     == Global variables == Line 32  C     == Global variables ==
32  #include "GMREDI.h"  #include "GMREDI.h"
33  #include "PARAMS.h"  #include "PARAMS.h"
34    
35    #ifdef ALLOW_AUTODIFF_TAMC
36    #include "tamc.h"
37    #include "tamc_keys.h"
38    #endif /* ALLOW_AUTODIFF_TAMC */
39    
40  C     == Routine arguments ==  C     == Routine arguments ==
41  C  C
42          _RL taperX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43          _RL taperY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46        _RL dSigmaDrW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDrW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47        _RL dSigmaDrS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDrS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48        _RL taperX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49        _RL taperY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50        _RL depthZ        _RS depthZ
51        INTEGER bi,bj,myThid        INTEGER K,bi,bj,myThid
52  CEndOfInterface  CEndOfInterface
53    
54  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
55    #ifdef GM_BOLUS_ADVEC
56    
57  C     == Local variables ==  C     == Local variables ==
58        _RL Small_Number        _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59        PARAMETER(Small_Number=1.D-12)        _RL f1,Smod,f2,Rnondim
       _RL gradSmod,f1,Smod,f2,Rnondim,Cspd,Lrho  
       _RL dSigmaDrLtd, dRdSigmaLtd  
60        _RL maxSlopeSqr        _RL maxSlopeSqr
61          _RL slopeCutoff
62        _RL fpi        _RL fpi
63        PARAMETER(fpi=3.141592653589793047592d0)        PARAMETER(fpi=3.141592653589793047592d0)
64        INTEGER i,j        INTEGER i,j
65    #ifdef GMREDI_WITH_STABLE_ADJOINT
66          _RL slopeTmpSpec,slopeMaxSpec
67    #endif
68    
69    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
70    
71          slopeCutoff = SQRT( GM_slopeSqCutoff )
72    
73    #ifdef ALLOW_AUTODIFF_TAMC
74          act1 = bi - myBxLo(myThid)
75          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
76          act2 = bj - myByLo(myThid)
77          max2 = myByHi(myThid) - myByLo(myThid) + 1
78          act3 = myThid - 1
79          max3 = nTx*nTy
80          act4 = ikey_dynamics - 1
81          igmkey = (act1 + 1) + act2*max1
82         &                    + act3*max1*max2
83         &                    + act4*max1*max2*max3
84          kkey = (igmkey-1)*Nr + k
85    #endif /* ALLOW_AUTODIFF_TAMC */
86    
87    #ifndef GMREDI_WITH_STABLE_ADJOINT
88    c common case:
89    
90        IF (GM_taper_scheme.EQ.'orig' .OR.        IF (GM_taper_scheme.EQ.'orig' .OR.
91       &    GM_taper_scheme.EQ.'clipping') THEN       &    GM_taper_scheme.EQ.'clipping') THEN
92    
93    #ifdef GM_EXCLUDE_CLIPPING
94    
95            STOP 'Need to compile without "#define GM_EXCLUDE_CLIPPING"'
96    
97    #else  /* GM_EXCLUDE_CLIPPING */
98    
99  C-      Original implementation in mitgcmuv  C-      Original implementation in mitgcmuv
100  C       (this turns out to be the same as Cox slope clipping)  C       (this turns out to be the same as Cox slope clipping)
101    
102  C-      Cox 1987 "Slope clipping"  C-- X-comp
         DO j=1-Oly+1,sNy+Oly-1  
          DO i=1-Olx+1,sNx+Olx-1  
103    
104  c         gradSmod=SlopeX(i,j)*SlopeX(i,j)  #ifdef ALLOW_AUTODIFF_TAMC
105  c    &            +SlopeY(i,j)*SlopeY(i,j)          DO j=1-Oly,sNy+Oly
106  c         if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)           DO i=1-Olx+1,sNx+Olx
107            gradSmod=abs(SlopeX(i,j))            dSigmaDrLtd(i,j) = 0. _d 0
108             ENDDO
109            dSigmaDrLtd = -(Small_Number+gradSmod*GM_rMaxSlope)          ENDDO
110            IF (dSigmaDrW(i,j).GE.dSigmaDrLtd)  #endif /* ALLOW_AUTODIFF_TAMC */
      &        dSigmaDrW(i,j) = dSigmaDrLtd  
           SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)  
111    
112            gradSmod=abs(SlopeY(i,j))  C-      Cox 1987 "Slope clipping"
113            dSigmaDrLtd = -(Small_Number+gradSmod*GM_rMaxSlope)          DO j=1-Oly,sNy+Oly
114            IF (dSigmaDrS(i,j).GE.dSigmaDrLtd)           DO i=1-Olx+1,sNx+Olx
115       &        dSigmaDrS(i,j) = dSigmaDrLtd            dSigmaDrLtd(i,j) = -(GM_Small_Number+
116            SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)       &     ABS(SlopeX(i,j))*GM_rMaxSlope)
117             ENDDO
118            ENDDO
119    #ifdef ALLOW_AUTODIFF_TAMC
120    CADJ STORE dSigmaDrLtd(:,:)  = comlev1_bibj_k, key=kkey, byte=isbyte
121    CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
122    #endif
123            DO j=1-Oly,sNy+Oly
124             DO i=1-Olx+1,sNx+Olx
125              IF (dSigmaDrW(i,j).GE.dSigmaDrLtd(i,j))
126         &        dSigmaDrW(i,j) = dSigmaDrLtd(i,j)
127             ENDDO
128            ENDDO
129    #ifdef ALLOW_AUTODIFF_TAMC
130    CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
131    #endif
132            DO j=1-Oly,sNy+Oly
133             DO i=1-Olx+1,sNx+Olx
134              SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
135              taperX(i,j) = 1. _d 0
136             ENDDO
137            ENDDO
138    
139            taperX(i,j)=1. _d 0  C-- Y-comp
           taperY(i,j)=1. _d 0  
140    
141    #ifdef ALLOW_AUTODIFF_TAMC
142            DO j=1-Oly+1,sNy+Oly
143             DO i=1-Olx,sNx+Olx
144              dSigmaDrLtd(i,j) = 0. _d 0
145             ENDDO
146            ENDDO
147    #endif /* ALLOW_AUTODIFF_TAMC */
148            DO j=1-Oly+1,sNy+Oly
149             DO i=1-Olx,sNx+Olx
150              dSigmaDrLtd(i,j) = -(GM_Small_Number+
151         &     ABS(SlopeY(i,j))*GM_rMaxSlope)
152           ENDDO           ENDDO
153          ENDDO          ENDDO
154    #ifdef ALLOW_AUTODIFF_TAMC
155    CADJ STORE dSigmaDrLtd(:,:)  = comlev1_bibj_k, key=kkey, byte=isbyte
156    CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
157    #endif
158            DO j=1-Oly+1,sNy+Oly
159             DO i=1-Olx,sNx+Olx
160              IF (dSigmaDrS(i,j).GE.dSigmaDrLtd(i,j))
161         &        dSigmaDrS(i,j) = dSigmaDrLtd(i,j)
162             ENDDO
163            ENDDO
164    #ifdef ALLOW_AUTODIFF_TAMC
165    CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
166    #endif
167            DO j=1-Oly+1,sNy+Oly
168             DO i=1-Olx,sNx+Olx
169              SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
170              taperY(i,j) = 1. _d 0
171             ENDDO
172            ENDDO
173    
174    #endif /* GM_EXCLUDE_CLIPPING */
175    
176          ELSEIF (GM_taper_scheme.EQ.'fm07') THEN
177    
178            STOP 'GMREDI_SLOPE_PSI: AdvForm not yet implemented for fm07'
179    
180        ELSE        ELSE
181    
182  C- Compute the slope, no clipping, but avoid reverse slope in negatively  #ifdef GM_EXCLUDE_TAPERING
 C                                  stratified (Sigma_Z > 0) region :  
         DO j=1-Oly+1,sNy+Oly-1  
          DO i=1-Olx+1,sNx+Olx-1  
183    
184            dSigmaDrLtd = -Small_Number          STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
           IF (dSigmaDrW(i,j).GE.dSigmaDrLtd)  
      &        dSigmaDrW(i,j) = dSigmaDrLtd  
           dRdSigmaLtd = 1./dSigmaDrW(i,j)  
           SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)  
185    
186            dSigmaDrLtd = -Small_Number  #else  /* GM_EXCLUDE_TAPERING */
           IF (dSigmaDrS(i,j).GE.dSigmaDrLtd)  
      &        dSigmaDrS(i,j) = dSigmaDrLtd  
           SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)  
187    
188  c         SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)  #ifdef ALLOW_AUTODIFF_TAMC
189  c    &                 +SlopeY(i,j)*SlopeY(i,j)  CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
190    CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
191    #endif
192    
193            taperX(i,j)=1. _d 0  C- Compute the slope, no clipping, but avoid reverse slope in negatively
194            taperY(i,j)=1. _d 0  C                                  stratified (Sigma_Z > 0) region :
195            DO j=1-Oly,sNy+Oly
196             DO i=1-Olx+1,sNx+Olx
197              IF (dSigmaDrW(i,j).GE.-GM_Small_Number)
198         &        dSigmaDrW(i,j) = -GM_Small_Number
199             ENDDO
200            ENDDO
201    #ifdef ALLOW_AUTODIFF_TAMC
202    CADJ STORE dsigmadrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
203    #endif
204            DO j=1-Oly,sNy+Oly
205             DO i=1-Olx+1,sNx+Olx
206              SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
207              taperX(i,j) = 1. _d 0
208             ENDDO
209            ENDDO
210    #ifdef ALLOW_AUTODIFF_TAMC
211    CADJ STORE slopex(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
212    #endif
213            DO j=1-Oly,sNy+Oly
214             DO i=1-Olx+1,sNx+Olx
215              IF ( ABS(SlopeX(i,j)) .GE. slopeCutoff ) THEN
216                 SlopeX(i,j) = SIGN(slopeCutoff,SlopeX(i,j))
217                 taperX(i,j) = 0. _d 0
218              ENDIF
219             ENDDO
220            ENDDO
221    
222    #ifdef ALLOW_AUTODIFF_TAMC
223    CADJ STORE slopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
224    CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
225    #endif
226    
227            DO j=1-Oly+1,sNy+Oly
228             DO i=1-Olx,sNx+Olx
229              IF (dSigmaDrS(i,j).GE.-GM_Small_Number)
230         &        dSigmaDrS(i,j) = -GM_Small_Number
231             ENDDO
232            ENDDO
233    #ifdef ALLOW_AUTODIFF_TAMC
234    CADJ STORE dsigmadrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
235    #endif
236            DO j=1-Oly+1,sNy+Oly
237             DO i=1-Olx,sNx+Olx
238              SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
239              taperY(i,j) = 1. _d 0
240             ENDDO
241            ENDDO
242    #ifdef ALLOW_AUTODIFF_TAMC
243    CADJ STORE slopey(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
244    #endif
245            DO j=1-Oly+1,sNy+Oly
246             DO i=1-Olx,sNx+Olx
247              IF ( ABS(SlopeY(i,j)) .GE. slopeCutoff ) THEN
248                 SlopeY(i,j) = SIGN(slopeCutoff,SlopeY(i,j))
249                 taperY(i,j) = 0. _d 0
250              ENDIF
251           ENDDO           ENDDO
252          ENDDO          ENDDO
253    
254  C- Compute the tapering function for the GM+Redi tensor :  C- Compute the tapering function for the GM+Redi tensor :
255    
256    #ifdef ALLOW_AUTODIFF_TAMC
257    CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
258    CADJ STORE slopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
259    #endif
260    
261         IF (GM_taper_scheme.EQ.'linear') THEN         IF (GM_taper_scheme.EQ.'linear') THEN
262    
263  C-      Simplest adiabatic tapering = Smax/Slope (linear)  C-      Simplest adiabatic tapering = Smax/Slope (linear)
264          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly,sNy+Oly
265           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx
266              Smod = ABS(SlopeX(i,j))
267            IF (abs(SlopeX(i,j)).GT.GM_maxSlope)            IF ( Smod .GT. GM_maxSlope .AND.
268       &      taperX(i,j)=GM_maxSlope/abs(SlopeX(i,j))       &           Smod .LT. slopeCutoff )
269            IF (abs(SlopeY(i,j)).GT.GM_maxSlope)       &           taperX(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
270       &      taperY(i,j)=GM_maxSlope/abs(SlopeY(i,j))           ENDDO
271            ENDDO
272            DO j=1-Oly+1,sNy+Oly
273             DO i=1-Olx,sNx+Olx
274              Smod = ABS(SlopeY(i,j))
275              IF ( Smod .GT. GM_maxSlope .AND.
276         &           Smod .LT. slopeCutoff )
277         &           taperY(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
278           ENDDO           ENDDO
279          ENDDO          ENDDO
280    
# Line 135  C-      Simplest adiabatic tapering = Sm Line 282  C-      Simplest adiabatic tapering = Sm
282    
283  C-      Gerdes, Koberle and Willebrand, Clim. Dyn. 1991  C-      Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
284          maxSlopeSqr = GM_maxSlope*GM_maxSlope          maxSlopeSqr = GM_maxSlope*GM_maxSlope
285          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly,sNy+Oly
286           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx
287              IF ( ABS(SlopeX(i,j)) .GT. GM_maxSlope .AND.
288            IF (abs(SlopeX(i,j)).GT.GM_maxSlope)       &           ABS(SlopeX(i,j)) .LT. slopeCutoff )
289       &      taperX(i,j)=maxSlopeSqr/(SlopeX(i,j)*SlopeX(i,j))       &           taperX(i,j)=maxSlopeSqr/
290            IF (abs(SlopeY(i,j)).GT.GM_maxSlope)       &           ( SlopeX(i,j)*SlopeX(i,j) + GM_Small_Number )
291       &      taperY(i,j)=maxSlopeSqr/(SlopeY(i,j)*SlopeY(i,j))           ENDDO
292            ENDDO
293            DO j=1-Oly+1,sNy+Oly
294             DO i=1-Olx,sNx+Olx
295              IF ( ABS(SlopeY(i,j)) .GT. GM_maxSlope .AND.
296         &           ABS(SlopeY(i,j)) .LT. slopeCutoff )
297         &           taperY(i,j)=maxSlopeSqr/
298         &           ( SlopeY(i,j)*SlopeY(i,j) + GM_Small_Number )
299           ENDDO           ENDDO
300          ENDDO          ENDDO
301    
302         ELSEIF (GM_taper_scheme.EQ.'dm95') THEN         ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
303    
304  C-      Danabasoglu and McWilliams, J. Clim. 1995  C-      Danabasoglu and McWilliams, J. Clim. 1995
305          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly,sNy+Oly
306           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx
307              Smod = ABS(SlopeX(i,j))
308            Smod = abs(SlopeX(i,j))            taperX(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
309            taperX(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))           ENDDO
310            Smod = abs(SlopeY(i,j))          ENDDO
311            taperY(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))          DO j=1-Oly+1,sNy+Oly
312             DO i=1-Olx,sNx+Olx
313              Smod = ABS(SlopeY(i,j))
314              taperY(i,j)=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
315           ENDDO           ENDDO
316          ENDDO          ENDDO
317    
318         ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN         ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
319    
320  C-      Large, Danabasoglu and Doney, JPO 1997  C-      Large, Danabasoglu and Doney, JPO 1997
         DO j=1-Oly+1,sNy+Oly-1  
          DO i=1-Olx+1,sNx+Olx-1  
321    
322            Cspd=2.          DO j=1-Oly,sNy+Oly
323            Lrho=100.e3           DO i=1-Olx+1,sNx+Olx
324            if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))            Smod = ABS(SlopeX(i,j))
325            Lrho=min(Lrho,100.e3)            IF ( Smod .LT. slopeCutoff ) THEN
326            Lrho=max(Lrho,15.e3)              f1=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
327                IF (Smod.NE.0.) THEN
328            Smod = abs(SlopeX(i,j))                Rnondim = -depthZ/(LrhoW(i,j)*Smod)
329            f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))              ELSE
330            if (Smod.NE.0.) then                Rnondim = 1.
331              Rnondim=depthZ/(Lrho*Smod)              ENDIF
332            else              IF ( Rnondim.GE.1. _d 0 ) THEN
333              Rnondim=0.                f2 = 1. _d 0
334            endif              ELSE
335            f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))                f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ))
336            taperX(i,j)=f1*f2              ENDIF
337                taperX(i,j)=f1*f2
338            Smod = abs(SlopeY(i,j))            ENDIF
339            f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))           ENDDO
340            if (Smod.NE.0.) then          ENDDO
             Rnondim=depthZ/(Lrho*Smod)  
           else  
             Rnondim=0.  
           endif  
           f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))  
           taperY(i,j)=f1*f2  
341    
342            DO j=1-Oly+1,sNy+Oly
343             DO i=1-Olx,sNx+Olx
344              Smod = ABS(SlopeY(i,j))
345              IF ( Smod .LT. slopeCutoff ) THEN
346                f1=op5*( 1. _d 0 + TANH( (GM_Scrit-Smod)/GM_Sd ))
347                IF (Smod.NE.0.) THEN
348                  Rnondim = -depthZ/(LrhoS(i,j)*Smod)
349                ELSE
350                  Rnondim = 1.
351                ENDIF
352                IF ( Rnondim.GE.1. _d 0 ) THEN
353                  f2 = 1. _d 0
354                ELSE
355                  f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ))
356                ENDIF
357                taperY(i,j)=f1*f2
358              ENDIF
359           ENDDO           ENDDO
360          ENDDO          ENDDO
361    
# Line 199  C-      Large, Danabasoglu and Doney, JP Line 363  C-      Large, Danabasoglu and Doney, JP
363          STOP 'GMREDI_SLOPE_PSI: Bad GM_taper_scheme'          STOP 'GMREDI_SLOPE_PSI: Bad GM_taper_scheme'
364         ENDIF         ENDIF
365    
366    #endif /* GM_EXCLUDE_TAPERING */
367    
368        ENDIF        ENDIF
369    
370    
371    #else  /* GMREDI_WITH_STABLE_ADJOINT */
372    c special choice for adjoint/optimization of parameters
373    c (~ strong clipping, reducing non linearity of psi=f(K))
374    
375    
376    #ifdef ALLOW_AUTODIFF_TAMC
377    CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
378    CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
379    #endif
380            DO j=1-Oly,sNy+Oly
381             DO i=1-Olx+1,sNx+Olx
382              IF (dSigmaDrW(i,j).GE.-GM_Small_Number)
383         &        dSigmaDrW(i,j) = -GM_Small_Number
384             ENDDO
385            ENDDO
386    #ifdef ALLOW_AUTODIFF_TAMC
387    CADJ STORE dsigmadrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
388    #endif
389            DO j=1-Oly,sNy+Oly
390             DO i=1-Olx+1,sNx+Olx
391              SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
392              taperX(i,j) = 1. _d 0
393             ENDDO
394            ENDDO
395    
396    #ifdef ALLOW_AUTODIFF_TAMC
397    CADJ STORE slopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
398    CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
399    #endif
400    
401            DO j=1-Oly+1,sNy+Oly
402             DO i=1-Olx,sNx+Olx
403              IF (dSigmaDrS(i,j).GE.-GM_Small_Number)
404         &        dSigmaDrS(i,j) = -GM_Small_Number
405             ENDDO
406            ENDDO
407    #ifdef ALLOW_AUTODIFF_TAMC
408    CADJ STORE dsigmadrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
409    #endif
410            DO j=1-Oly+1,sNy+Oly
411             DO i=1-Olx,sNx+Olx
412              SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
413              taperY(i,j) = 1. _d 0
414             ENDDO
415            ENDDO
416    
417            slopeMaxSpec=1. _d -4
418    
419    CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
420    CADJ STORE slopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
421    
422            DO j=1-Oly,sNy+Oly
423             DO i=1-Olx+1,sNx+Olx
424           slopeTmpSpec=ABS(SlopeX(i,j))
425           IF ( slopeTmpSpec .GT. slopeMaxSpec ) then
426            SlopeX(i,j)=5.*SlopeX(i,j)*slopeMaxSpec/slopeTmpSpec
427           ELSE
428            SlopeX(i,j)=5.*SlopeX(i,j)
429           ENDIF
430           taperX(i,j)=1.
431             ENDDO
432            ENDDO
433            DO j=1-Oly+1,sNy+Oly
434             DO i=1-Olx,sNx+Olx
435           slopeTmpSpec=ABS(SlopeY(i,j))
436           IF ( slopeTmpSpec .GT. slopeMaxSpec ) then
437            SlopeY(i,j)=5.*SlopeY(i,j)*slopeMaxSpec/slopeTmpSpec
438           ELSE
439            SlopeY(i,j)=5.*SlopeY(i,j)
440           ENDIF
441           taperY(i,j)=1.
442             ENDDO
443            ENDDO
444    
445    #endif /* GMREDI_WITH_STABLE_ADJOINT */
446    
447    #endif /* BOLUS_ADVEC */
448  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
449    
450        RETURN        RETURN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22