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