/[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.6 by heimbach, Tue Jan 21 19:34:13 2003 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(
8       I             dSigmaDrW,dSigmaDrS,       I             dSigmaDrW,dSigmaDrS,
9       I             depthZ,       I             depthZ,K,
10       U             SlopeX, SlopeY,       U             SlopeX, SlopeY,
11       O             taperX, taperY,       O             taperX, taperY,
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 conatins 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       conatins the height (m) of level |
21  C     | On exit:                                                 |  C     | On exit:                                                 |
22  C     |            dSigmaDrW conatins the effective dSig/dz   |  C     |            dSigmaDrW conatins 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 SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
# Line 41  C Line 46  C
46        _RL taperX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47        _RL taperY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48        _RL depthZ        _RL depthZ
49        INTEGER bi,bj,myThid        INTEGER K,bi,bj,myThid
50  CEndOfInterface  CEndOfInterface
51    
52  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
53    #ifdef GM_BOLUS_ADVEC
54    
55  C     == Local variables ==  C     == Local variables ==
56        _RL Small_Number        _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57        PARAMETER(Small_Number=1.D-12)        _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58        _RL gradSmod,f1,Smod,f2,Rnondim,Cspd,Lrho        _RL drdsigmaltd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59        _RL dSigmaDrLtd, dRdSigmaLtd        _RL  SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60          _RL f1,Smod,f2,Rnondim,Cspd,Lrho
61        _RL maxSlopeSqr        _RL maxSlopeSqr
62          _RL tmpvar
63        _RL fpi        _RL fpi
64        PARAMETER(fpi=3.141592653589793047592d0)        PARAMETER(fpi=3.141592653589793047592d0)
65        INTEGER i,j        INTEGER i,j
66    
67    #ifdef ALLOW_AUTODIFF_TAMC
68          act1 = bi - myBxLo(myThid)
69          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
70          act2 = bj - myByLo(myThid)
71          max2 = myByHi(myThid) - myByLo(myThid) + 1
72          act3 = myThid - 1
73          max3 = nTx*nTy
74          act4 = ikey_dynamics - 1
75          igmkey = (act1 + 1) + act2*max1
76         &                    + act3*max1*max2
77         &                    + act4*max1*max2*max3
78          kkey = (igmkey-1)*Nr + k
79    #endif /* ALLOW_AUTODIFF_TAMC */
80    
81        IF (GM_taper_scheme.EQ.'orig' .OR.        IF (GM_taper_scheme.EQ.'orig' .OR.
82       &    GM_taper_scheme.EQ.'clipping') THEN       &    GM_taper_scheme.EQ.'clipping') THEN
83    
84    #ifdef GM_EXCLUDE_CLIPPING
85    
86            STOP 'Need to compile without "#define GM_EXCLUDE_CLIPPING"'
87    
88    #else  /* GM_EXCLUDE_CLIPPING */
89    
90  C-      Original implementation in mitgcmuv  C-      Original implementation in mitgcmuv
91  C       (this turns out to be the same as Cox slope clipping)  C       (this turns out to be the same as Cox slope clipping)
92    
93    C-- X-comp
94    
95    #ifdef ALLOW_AUTODIFF_TAMC
96          DO j=1-Oly+1,sNy+Oly-1
97           DO i=1-Olx+1,sNx+Olx-1
98            dSigmaDrLtd(i,j) = 0. _d 0
99           ENDDO
100          ENDDO
101    #endif /* ALLOW_AUTODIFF_TAMC */
102    
103  C-      Cox 1987 "Slope clipping"  C-      Cox 1987 "Slope clipping"
104          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
105           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
106              dsigmadrltd(i,j) = -(GM_Small_Number+
107  c         gradSmod=SlopeX(i,j)*SlopeX(i,j)       &     abs(SlopeX(i,j))*GM_rMaxSlope)
108  c    &            +SlopeY(i,j)*SlopeY(i,j)           ENDDO
109  c         if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)          ENDDO
110            gradSmod=abs(SlopeX(i,j))  #ifdef ALLOW_AUTODIFF_TAMC
111    CADJ STORE dSigmaDrltd(:,:)  = comlev1_bibj_k, key=kkey, byte=isbyte
112            dSigmaDrLtd = -(Small_Number+gradSmod*GM_rMaxSlope)  CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
113            IF (dSigmaDrW(i,j).GE.dSigmaDrLtd)  #endif
114       &        dSigmaDrW(i,j) = dSigmaDrLtd          DO j=1-Oly+1,sNy+Oly-1
115             DO i=1-Olx+1,sNx+Olx-1
116              IF (dSigmaDrW(i,j).GE.dsigmadrltd(i,j))
117         &        dSigmaDrW(i,j) = dsigmadrltd(i,j)
118             ENDDO
119            ENDDO
120    #ifdef ALLOW_AUTODIFF_TAMC
121    CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
122    #endif
123            DO j=1-Oly+1,sNy+Oly-1
124             DO i=1-Olx+1,sNx+Olx-1
125            SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)            SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
126              taperX(i,j)=1. _d 0
127             ENDDO
128            ENDDO
129    
130            gradSmod=abs(SlopeY(i,j))  C-- Y-comp
           dSigmaDrLtd = -(Small_Number+gradSmod*GM_rMaxSlope)  
           IF (dSigmaDrS(i,j).GE.dSigmaDrLtd)  
      &        dSigmaDrS(i,j) = dSigmaDrLtd  
           SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)  
131    
132            taperX(i,j)=1. _d 0  #ifdef ALLOW_AUTODIFF_TAMC
133          DO j=1-Oly+1,sNy+Oly-1
134           DO i=1-Olx+1,sNx+Olx-1
135            dSigmaDrLtd(i,j) = 0. _d 0
136           ENDDO
137          ENDDO
138    #endif /* ALLOW_AUTODIFF_TAMC */
139            DO j=1-Oly+1,sNy+Oly-1
140             DO i=1-Olx+1,sNx+Olx-1
141              dsigmadrltd(i,j) = -(GM_Small_Number+
142         &     abs(SlopeY(i,j))*GM_rMaxSlope)
143             ENDDO
144            ENDDO
145    #ifdef ALLOW_AUTODIFF_TAMC
146    CADJ STORE dSigmaDrltd(:,:)  = comlev1_bibj_k, key=kkey, byte=isbyte
147    CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
148    #endif
149            DO j=1-Oly+1,sNy+Oly-1
150             DO i=1-Olx+1,sNx+Olx-1
151              IF (dSigmaDrS(i,j).GE.dsigmadrltd(i,j))
152         &        dSigmaDrS(i,j) = dsigmadrltd(i,j)
153             ENDDO
154            ENDDO
155    #ifdef ALLOW_AUTODIFF_TAMC
156    CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
157    #endif
158            DO j=1-Oly+1,sNy+Oly-1
159             DO i=1-Olx+1,sNx+Olx-1
160              SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
161            taperY(i,j)=1. _d 0            taperY(i,j)=1. _d 0
   
162           ENDDO           ENDDO
163          ENDDO          ENDDO
164    
165    #endif /* GM_EXCLUDE_CLIPPING */
166    
167        ELSE        ELSE
168    
169    #ifdef GM_EXCLUDE_TAPERING
170    
171            STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
172    
173    #else  /* GM_EXCLUDE_TAPERING */
174    
175    #ifdef ALLOW_AUTODIFF_TAMC
176    CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
177    CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
178    #endif
179    
180            tmpvar = sqrt( GM_slopeSqCutoff )
181    
182  C- Compute the slope, no clipping, but avoid reverse slope in negatively  C- Compute the slope, no clipping, but avoid reverse slope in negatively
183  C                                  stratified (Sigma_Z > 0) region :  C                                  stratified (Sigma_Z > 0) region :
184          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
185           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
186              dsigmadrltd(i,j) = -GM_Small_Number
187            dSigmaDrLtd = -Small_Number            IF (dSigmaDrW(i,j).GE.dsigmadrltd(i,j))
188            IF (dSigmaDrW(i,j).GE.dSigmaDrLtd)       &        dSigmaDrW(i,j) = dsigmadrltd(i,j)
189       &        dSigmaDrW(i,j) = dSigmaDrLtd           ENDDO
190            dRdSigmaLtd = 1./dSigmaDrW(i,j)          ENDDO
191    #ifdef ALLOW_AUTODIFF_TAMC
192    CADJ STORE dsigmadrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
193    #endif
194            DO j=1-Oly+1,sNy+Oly-1
195             DO i=1-Olx+1,sNx+Olx-1
196              drdsigmaltd(i,j) = 1./dSigmaDrW(i,j)
197            SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)            SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
198              taperX(i,j)= 1. _d 0
199             ENDDO
200            ENDDO
201    #ifdef ALLOW_AUTODIFF_TAMC
202    CADJ STORE slopex(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
203    #endif
204            DO j=1-Oly+1,sNy+Oly-1
205             DO i=1-Olx+1,sNx+Olx-1
206              IF ( ABS(SlopeX(i,j)) .GE. tmpvar ) THEN
207                 SlopeX(i,j) = SIGN(tmpvar,SlopeX(i,j))
208                 taperX(i,j) = 0. _d 0
209              ENDIF
210             ENDDO
211            ENDDO
212    
213            dSigmaDrLtd = -Small_Number  #ifdef ALLOW_AUTODIFF_TAMC
214            IF (dSigmaDrS(i,j).GE.dSigmaDrLtd)  CADJ STORE slopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
215       &        dSigmaDrS(i,j) = dSigmaDrLtd  CADJ STORE dSigmaDrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
216            SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)  #endif
   
 c         SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)  
 c    &                 +SlopeY(i,j)*SlopeY(i,j)  
217    
218            taperX(i,j)=1. _d 0          DO j=1-Oly+1,sNy+Oly-1
219             DO i=1-Olx+1,sNx+Olx-1
220              dsigmadrltd(i,j) = -GM_Small_Number
221              IF (dSigmaDrS(i,j).GE.dsigmadrltd(i,j))
222         &        dSigmaDrS(i,j) = dsigmadrltd(i,j)
223             ENDDO
224            ENDDO
225    #ifdef ALLOW_AUTODIFF_TAMC
226    CADJ STORE dsigmadrS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
227    #endif
228            DO j=1-Oly+1,sNy+Oly-1
229             DO i=1-Olx+1,sNx+Olx-1
230              SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
231            taperY(i,j)=1. _d 0            taperY(i,j)=1. _d 0
232             ENDDO
233            ENDDO
234    #ifdef ALLOW_AUTODIFF_TAMC
235    CADJ STORE slopey(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
236    #endif
237            DO j=1-Oly+1,sNy+Oly-1
238             DO i=1-Olx+1,sNx+Olx-1
239              IF ( ABS(SlopeY(i,j)) .GE. tmpvar ) THEN
240                 SlopeY(i,j) = SIGN(tmpvar,SlopeY(i,j))
241                 taperY(i,j) = 0. _d 0
242              ENDIF
243           ENDDO           ENDDO
244          ENDDO          ENDDO
245    
246  C- Compute the tapering function for the GM+Redi tensor :  C- Compute the tapering function for the GM+Redi tensor :
247    
248    #ifdef ALLOW_AUTODIFF_TAMC
249    CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
250    CADJ STORE slopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
251    #endif
252    
253         IF (GM_taper_scheme.EQ.'linear') THEN         IF (GM_taper_scheme.EQ.'linear') THEN
254    
255  C-      Simplest adiabatic tapering = Smax/Slope (linear)  C-      Simplest adiabatic tapering = Smax/Slope (linear)
256          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
257           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
258    
259            IF (abs(SlopeX(i,j)).GT.GM_maxSlope)            Smod = abs(SlopeX(i,j))
260       &      taperX(i,j)=GM_maxSlope/abs(SlopeX(i,j))            IF ( Smod .GT. GM_maxSlope .AND.
261            IF (abs(SlopeY(i,j)).GT.GM_maxSlope)       &           Smod .LT. tmpvar )
262       &      taperY(i,j)=GM_maxSlope/abs(SlopeY(i,j))       &           taperX(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
263              Smod = abs(SlopeY(i,j))
264              IF ( Smod .GT. GM_maxSlope .AND.
265         &           Smod .LT. tmpvar )
266         &           taperY(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
267    
268           ENDDO           ENDDO
269          ENDDO          ENDDO
# Line 134  C-      Simplest adiabatic tapering = Sm Line 271  C-      Simplest adiabatic tapering = Sm
271         ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN         ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
272    
273  C-      Gerdes, Koberle and Willebrand, Clim. Dyn. 1991  C-      Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
274    
275          maxSlopeSqr = GM_maxSlope*GM_maxSlope          maxSlopeSqr = GM_maxSlope*GM_maxSlope
276          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
277           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
278    
279            IF (abs(SlopeX(i,j)).GT.GM_maxSlope)            IF ( abs(SlopeX(i,j)) .GT. GM_maxSlope .AND.
280       &      taperX(i,j)=maxSlopeSqr/(SlopeX(i,j)*SlopeX(i,j))       &           abs(SlopeX(i,j)) .LT. tmpvar )
281            IF (abs(SlopeY(i,j)).GT.GM_maxSlope)       &           taperX(i,j)=maxSlopeSqr/
282       &      taperY(i,j)=maxSlopeSqr/(SlopeY(i,j)*SlopeY(i,j))       &           ( SlopeX(i,j)*SlopeX(i,j) + GM_Small_Number )
283              IF ( abs(SlopeY(i,j)) .GT. GM_maxSlope .AND.
284         &           abs(SlopeY(i,j)) .LT. tmpvar )
285         &           taperY(i,j)=maxSlopeSqr/
286         &           ( SlopeY(i,j)*SlopeY(i,j) + GM_Small_Number )
287    
288           ENDDO           ENDDO
289          ENDDO          ENDDO
# Line 153  C-      Danabasoglu and McWilliams, J. C Line 295  C-      Danabasoglu and McWilliams, J. C
295           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
296    
297            Smod = abs(SlopeX(i,j))            Smod = abs(SlopeX(i,j))
298            taperX(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))            taperX(i,j)=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
299            Smod = abs(SlopeY(i,j))            Smod = abs(SlopeY(i,j))
300            taperY(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))            taperY(i,j)=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
301    
302           ENDDO           ENDDO
303          ENDDO          ENDDO
# Line 166  C-      Large, Danabasoglu and Doney, JP Line 308  C-      Large, Danabasoglu and Doney, JP
308          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
309           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
310    
311            Cspd=2.            Cspd=2. _d 0
312            Lrho=100.e3            Lrho=100. _d 3
313            if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))            if (fCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(fCori(i,j,bi,bj))
314            Lrho=min(Lrho,100.e3)            Lrho=min(Lrho , 100. _d 3)
315            Lrho=max(Lrho,15.e3)            Lrho=max(Lrho , 15. _d 3)
316    
317            Smod = abs(SlopeX(i,j))            Smod = abs(SlopeX(i,j))
318            f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))            if ( Smod .LT. tmpvar ) then
319              f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
320            if (Smod.NE.0.) then            if (Smod.NE.0.) then
321              Rnondim=depthZ/(Lrho*Smod)              Rnondim=depthZ/(Lrho*Smod)
322            else            else
323              Rnondim=0.              Rnondim=0.
324            endif            endif
325            f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))            f2=op5*( 1. _d 0 + sin( fpi*(Rnondim-op5)))
326            taperX(i,j)=f1*f2            taperX(i,j)=f1*f2
327              endif
328    
329            Smod = abs(SlopeY(i,j))            Smod = abs(SlopeY(i,j))
330            f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))            if ( Smod .LT. tmpvar ) then
331              f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
332            if (Smod.NE.0.) then            if (Smod.NE.0.) then
333              Rnondim=depthZ/(Lrho*Smod)              Rnondim=depthZ/(Lrho*Smod)
334            else            else
335              Rnondim=0.              Rnondim=0.
336            endif            endif
337            f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))            f2=op5*( 1. _d 0 + sin( fpi*(Rnondim-op5)))
338            taperY(i,j)=f1*f2            taperY(i,j)=f1*f2
339              endif
340    
341           ENDDO           ENDDO
342          ENDDO          ENDDO
# Line 199  C-      Large, Danabasoglu and Doney, JP Line 345  C-      Large, Danabasoglu and Doney, JP
345          STOP 'GMREDI_SLOPE_PSI: Bad GM_taper_scheme'          STOP 'GMREDI_SLOPE_PSI: Bad GM_taper_scheme'
346         ENDIF         ENDIF
347    
348    #endif /* GM_EXCLUDE_TAPERING */
349    
350        ENDIF        ENDIF
351    
352    
353    #endif /* BOLUS_ADVEC */
354  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
355    
356        RETURN        RETURN
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22