/[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.4 by heimbach, Fri Jan 10 23:41:15 2003 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,K,  
      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 39  C     == Global variables == Line 39  C     == Global variables ==
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          _RS depthZ
52        INTEGER K,bi,bj,myThid        INTEGER K,bi,bj,myThid
53  CEndOfInterface  CEndOfInterface
54    
# Line 53  CEndOfInterface Line 56  CEndOfInterface
56  #ifdef GM_BOLUS_ADVEC  #ifdef GM_BOLUS_ADVEC
57    
58  C     == Local variables ==  C     == Local variables ==
       _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)  
59        _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60        _RL drdsigmaltd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL f1,Smod,f2,Rnondim
       _RL f1,Smod,f2,Rnondim,Cspd,Lrho  
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  #ifdef ALLOW_AUTODIFF_TAMC
70        act1 = bi - myBxLo(myThid)        act1 = bi - myBxLo(myThid)
71        max1 = myBxHi(myThid) - myBxLo(myThid) + 1        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
# Line 70  C     == Local variables == Line 74  C     == Local variables ==
74        act3 = myThid - 1        act3 = myThid - 1
75        max3 = nTx*nTy        max3 = nTx*nTy
76        act4 = ikey_dynamics - 1        act4 = ikey_dynamics - 1
77        ikey = (act1 + 1) + act2*max1        igmkey = (act1 + 1) + act2*max1
78       &                  + act3*max1*max2       &                    + act3*max1*max2
79       &                  + act4*max1*max2*max3       &                    + act4*max1*max2*max3
80        kkey = (ikey-1)*Nr + k        kkey = (igmkey-1)*Nr + k
81  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
82    
83  #ifdef ALLOW_AUTODIFF_TAMC  #ifndef GMREDI_WITH_STABLE_ADJOINT
84  CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  c common case:
 CADJ STORE slopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif  
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-- X-comp  C-- X-comp
99    
100  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
101        DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly,sNy+Oly
102         DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx
103          dSigmaDrLtd(i,j) = 0. _d 0            dSigmaDrLtd(i,j) = 0. _d 0
104         ENDDO           ENDDO
105        ENDDO          ENDDO
106  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
107    
108  C-      Cox 1987 "Slope clipping"  C-      Cox 1987 "Slope clipping"
109          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly,sNy+Oly
110           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx
111            dsigmadrltd(i,j) = -(GM_Small_Number+            dSigmaDrLtd(i,j) = -(GM_Small_Number+
112       &     abs(SlopeX(i,j))*GM_rMaxSlope)       &     ABS(SlopeX(i,j))*GM_rMaxSlope)
113           ENDDO           ENDDO
114          ENDDO          ENDDO
115  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
116  CADJ STORE dSigmaDrltd(:,:)  = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDrLtd(:,:)  = comlev1_bibj_k, key=kkey, byte=isbyte
117  CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
118  #endif  #endif
119          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly,sNy+Oly
120           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx
121            IF (dSigmaDrW(i,j).GE.dsigmadrltd(i,j))            IF (dSigmaDrW(i,j).GE.dSigmaDrLtd(i,j))
122       &        dSigmaDrW(i,j) = dsigmadrltd(i,j)       &        dSigmaDrW(i,j) = dSigmaDrLtd(i,j)
123           ENDDO           ENDDO
124          ENDDO          ENDDO
125  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
126  CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
127  #endif  #endif
128          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly,sNy+Oly
129           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx
130            SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)            SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
131            taperX(i,j)=1. _d 0            taperX(i,j) = 1. _d 0
132           ENDDO           ENDDO
133          ENDDO          ENDDO
134    
135  C-- Y-comp  C-- Y-comp
136    
137  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
138        DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly
139         DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx,sNx+Olx
140          dSigmaDrLtd(i,j) = 0. _d 0            dSigmaDrLtd(i,j) = 0. _d 0
141         ENDDO           ENDDO
142        ENDDO          ENDDO
143  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
144          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly
145           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx,sNx+Olx
146            dsigmadrltd(i,j) = -(GM_Small_Number+            dSigmaDrLtd(i,j) = -(GM_Small_Number+
147       &     abs(SlopeY(i,j))*GM_rMaxSlope)       &     ABS(SlopeY(i,j))*GM_rMaxSlope)
148           ENDDO           ENDDO
149          ENDDO          ENDDO
150  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
151  CADJ STORE dSigmaDrltd(:,:)  = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDrLtd(:,:)  = comlev1_bibj_k, key=kkey, byte=isbyte
152  CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
153  #endif  #endif
154          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly
155           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx,sNx+Olx
156            IF (dSigmaDrS(i,j).GE.dsigmadrltd(i,j))            IF (dSigmaDrS(i,j).GE.dSigmaDrLtd(i,j))
157       &        dSigmaDrS(i,j) = dsigmadrltd(i,j)       &        dSigmaDrS(i,j) = dSigmaDrLtd(i,j)
158           ENDDO           ENDDO
159          ENDDO          ENDDO
160  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
161  CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
162  #endif  #endif
163          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly
164           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx,sNx+Olx
165            SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)            SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
166            taperY(i,j)=1. _d 0            taperY(i,j) = 1. _d 0
167           ENDDO           ENDDO
168          ENDDO          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    #ifdef GM_EXCLUDE_TAPERING
179    
180            STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
181    
182    #else  /* GM_EXCLUDE_TAPERING */
183    
184  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
185  CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
186  CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 168  CADJ STORE dSigmaDrW(:,:)    = comlev1_b Line 188  CADJ STORE dSigmaDrW(:,:)    = comlev1_b
188    
189  C- Compute the slope, no clipping, but avoid reverse slope in negatively  C- Compute the slope, no clipping, but avoid reverse slope in negatively
190  C                                  stratified (Sigma_Z > 0) region :  C                                  stratified (Sigma_Z > 0) region :
191          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly,sNy+Oly
192           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx
193            dsigmadrltd(i,j) = -GM_Small_Number            IF (dSigmaDrW(i,j).GE.-GM_Small_Number)
194            IF (dSigmaDrW(i,j).GE.dsigmadrltd(i,j))       &        dSigmaDrW(i,j) = -GM_Small_Number
      &        dSigmaDrW(i,j) = dsigmadrltd(i,j)  
195           ENDDO           ENDDO
196          ENDDO          ENDDO
197  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
198  CADJ STORE dsigmadrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dsigmadrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
199  #endif  #endif
200          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly,sNy+Oly
201           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx
           drdsigmaltd(i,j) = 1./dSigmaDrW(i,j)  
202            SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)            SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
203            taperX(i,j)=1. _d 0            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           ENDDO
216          ENDDO          ENDDO
217    
# Line 191  CADJ STORE slopeY(:,:)       = comlev1_b Line 220  CADJ STORE slopeY(:,:)       = comlev1_b
220  CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
221  #endif  #endif
222    
223          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly
224           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx,sNx+Olx
225            dsigmadrltd(i,j) = -GM_Small_Number            IF (dSigmaDrS(i,j).GE.-GM_Small_Number)
226            IF (dSigmaDrS(i,j).GE.dsigmadrltd(i,j))       &        dSigmaDrS(i,j) = -GM_Small_Number
      &        dSigmaDrS(i,j) = dsigmadrltd(i,j)  
227           ENDDO           ENDDO
228          ENDDO          ENDDO
229  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
230  CADJ STORE dsigmadrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dsigmadrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
231  #endif  #endif
232          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly
233           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx,sNx+Olx
234            SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)            SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
235            taperY(i,j)=1. _d 0            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    
# Line 218  CADJ STORE slopeY(:,:)       = comlev1_b Line 257  CADJ STORE slopeY(:,:)       = comlev1_b
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 233  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)=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))           ENDDO
306            Smod = abs(SlopeY(i,j))          ENDDO
307            taperY(i,j)=op5*( 1. _d 0 + 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. _d 0          DO j=1-Oly,sNy+Oly
319            Lrho=100. _d 0           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=op5*( 1. _d 0 + 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=op5*( 1. _d 0 + sin( fpi*(Rnondim-op5)))                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=op5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))           ENDDO
336            if (Smod.NE.0.) then          ENDDO
             Rnondim=depthZ/(Lrho*Smod)  
           else  
             Rnondim=0.  
           endif  
           f2=op5*(1.+sin( fpi*(Rnondim-op5)))  
           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 297  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 */  #endif /* BOLUS_ADVEC */
444  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
445    
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22