/[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.4 by adcroft, Tue May 29 14:01:37 2001 UTC revision 1.12 by heimbach, Sun Mar 24 02:33:16 2002 UTC
# Line 1  Line 1 
 C $Header$  
 C $Name$  
1    
2  #include "GMREDI_OPTIONS.h"  #include "GMREDI_OPTIONS.h"
3    
4  CStartOfInterface  CStartOfInterface
5        SUBROUTINE GMREDI_SLOPE_LIMIT(        SUBROUTINE GMREDI_SLOPE_LIMIT(
6       I             dSigmadRReal,       I             dSigmaDrReal,
7       I             depthZ,       I             depthZ,
8       U             SlopeX, SlopeY,       U             SlopeX, SlopeY,
9       O             dRdSigmaLtd,       O             SlopeSqr, taperFct,
10       I             bi,bj, myThid )       I             bi,bj, myThid )
11  C     /==========================================================\  C     /==========================================================\
12  C     | SUBROUTINE GMREDI_SLOPE_LIMIT                            |  C     | SUBROUTINE GMREDI_SLOPE_LIMIT                            |
13  C     | o Calculate slopes for use in GM/Redi tensor             |  C     | o Calculate slopes for use in GM/Redi tensor             |
14  C     |==========================================================|  C     |==========================================================|
15  C     | On entry:                                                |  C     | On entry:                                                |
16  C     |            dSigmadRReal conatins the d/dz Sigma          |  C     |            dSigmaDrReal conatins the d/dz Sigma          |
17  C     |            SlopeX/Y     contains X/Y gradients of sigma  |  C     |            SlopeX/Y     contains X/Y gradients of sigma  |
18  C     |            depthZ       conatins the height (m) of level |  C     |            depthZ       conatins the height (m) of level |
19  C     | On exit:                                                 |  C     | On exit:                                                 |
20    C     |            dSigmaDrReal conatins the effective dSig/dz   |
21  C     |            SlopeX/Y     contains X/Y slopes              |  C     |            SlopeX/Y     contains X/Y slopes              |
22  C     |            dRdSigmaLtd  conatins the effective dSig/dz   |  C     |            SlopeSqr     contains Sx^2+Sy^2               |
23    C     |            taperFct     contains tapering funct. value ; |
24    C     |                         = 1 when using no tapering       |
25  C     \==========================================================/  C     \==========================================================/
26        IMPLICIT NONE        IMPLICIT NONE
27    
28  C     == Global variables ==  C     == Global variables ==
29  #include "SIZE.h"  #include "SIZE.h"
30    #include "GRID.h"
31  #include "EEPARAMS.h"  #include "EEPARAMS.h"
32  #include "GMREDI.h"  #include "GMREDI.h"
33  #include "PARAMS.h"  #include "PARAMS.h"
34    #ifdef ALLOW_AUTODIFF_TAMC
35    #include "tamc.h"
36    #include "tamc_keys.h"
37    #endif /* ALLOW_AUTODIFF_TAMC */
38    
39  C     == Routine arguments ==  C     == Routine arguments ==
40  C  C
41        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
42        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44        _RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45          _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46        _RL depthZ        _RL depthZ
47        INTEGER bi,bj,myThid        INTEGER bi,bj,myThid
48  CEndOfInterface  CEndOfInterface
# Line 44  CEndOfInterface Line 51  CEndOfInterface
51    
52  C     == Local variables ==  C     == Local variables ==
53        _RL Small_Number        _RL Small_Number
54        PARAMETER(Small_Number=1.E-12)        _RL Small_Taper
55        _RL gradSmod,stratLimit,f1,Smod,f2,Rnondim,Cspd,Lrho        PARAMETER(Small_Number=1.D-12)
56        _RL dSigmaDrLtd        PARAMETER(Small_Taper=1.D+03)
57          _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58          _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59          _RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60          _RL f1,Smod,f2,Rnondim,Cspd,Lrho
61          _RL maxSlopeSqr
62        _RL fpi        _RL fpi
63        PARAMETER(fpi=3.141592653589793047592)        PARAMETER(fpi=3.141592653589793047592d0)
64        INTEGER i,j        INTEGER i,j
       _RL     sx, sy  
65    
66        if (GM_taper_scheme.EQ.'orig') then  #ifdef ALLOW_AUTODIFF_TAMC
67          act1 = bi - myBxLo(myThid)
68          max1 = myBxHi(myThid) - myBxLo(myThid) + 1
69          act2 = bj - myByLo(myThid)
70          max2 = myByHi(myThid) - myByLo(myThid) + 1
71          act3 = myThid - 1
72          max3 = nTx*nTy
73          act4 = ikey_dynamics - 1
74          ikey = (act1 + 1) + act2*max1
75         &                  + act3*max1*max2
76         &                  + act4*max1*max2*max3
77    #endif /* ALLOW_AUTODIFF_TAMC */
78    
79          DO j=1-Oly+1,sNy+Oly-1
80           DO i=1-Olx+1,sNx+Olx-1
81            gradSmod(i,j)    = 0. _d 0
82            dSigmaDrLtd(i,j) = 0. _d 0
83           ENDDO
84          ENDDO
85    
86          IF (GM_taper_scheme.EQ.'orig' .OR.
87         &    GM_taper_scheme.EQ.'clipping') THEN
88    
89  C-      Original implementation in mitgcmuv  C-      Original implementation in mitgcmuv
90  C       (this turns out to be the same as Cox slope clipping)  C       (this turns out to be the same as Cox slope clipping)
91    
92    C-      Cox 1987 "Slope clipping"
93          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
94           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
95              IF ( SlopeX(i,j)*SlopeX(i,j) + SlopeY(i,j)*SlopeY(i,j)
96         &        .EQ. 0. ) THEN
97               gradSmod(i,j) = 0. _d 0
98              ELSE
99               gradSmod(i,j) = sqrt(SlopeX(i,j)*SlopeX(i,j)
100         &                         +SlopeY(i,j)*SlopeY(i,j))
101              ENDIF
102             ENDDO
103            ENDDO
104    
105    #ifdef ALLOW_AUTODIFF_TAMC
106    CADJ STORE gradSmod(:,:)     = comlev1_bibj, key=ikey, byte=isbyte
107    CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
108    #endif
109    
110            gradSmod=SlopeX(i,j)*SlopeX(i,j)          DO j=1-Oly+1,sNy+Oly-1
111       &            +SlopeY(i,j)*SlopeY(i,j)           DO i=1-Olx+1,sNx+Olx-1
112            if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)            IF (gradSmod(i,j) .NE. 0.) THEN
113               dSigmaDrLtd(i,j) = -gradSmod(i,j)*GM_rMaxSlope
114               IF (dSigmaDrReal(i,j) .GE. dSigmaDrLtd(i,j))
115         &         dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
116              ENDIF
117             ENDDO
118            ENDDO
119    
120            stratLimit=-Small_Number-gradSmod*GM_rMaxSlope  #ifdef ALLOW_AUTODIFF_TAMC
121            if (dSigmaDrReal(i,j).LT.stratLimit) then  CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
122              dSigmaDrLtd=dSigmaDrReal(i,j)  CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
123            else  CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
124              dSigmaDrLtd=stratLimit  #endif
           endif  
           dRdSigmaLtd(i,j)=1./dSigmaDrLtd  
           SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)  
           SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)  
125    
126            DO j=1-Oly+1,sNy+Oly-1
127             DO i=1-Olx+1,sNx+Olx-1
128              IF (gradSmod(i,j) .EQ. 0.) THEN
129               SlopeX(i,j) = 0. _d 0
130               SlopeY(i,j) = 0. _d 0
131              ELSE
132               dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)
133               SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)
134               SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)
135              ENDIF
136           ENDDO           ENDDO
137          ENDDO          ENDDO
138    
139        elseif (GM_taper_scheme.EQ.'clipping') then  #ifdef ALLOW_AUTODIFF_TAMC
140    CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
141    CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
142    #endif
143    
 C-      Cox 1987 "Slope clipping"  
144          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
145           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
146              SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)
147         &                 +SlopeY(i,j)*SlopeY(i,j)
148              taperFct(i,j)=1. _d 0
149             ENDDO
150            ENDDO
151    
152          ELSE
153    
154            gradSmod=SlopeX(i,j)*SlopeX(i,j)  C----------------------------------------------------------------------
      &            +SlopeY(i,j)*SlopeY(i,j)  
           if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)  
155    
156            dSigmaDrLtd=-(Small_Number+gradSmod*GM_rMaxSlope)  C- Compute the slope, no clipping, but avoid reverse slope in negatively
157            if (dSigmaDrReal(i,j).LT.dSigmaDrLtd)  C                                  stratified (Sigma_Z > 0) region :
      &            dSigmaDrLtd=dSigmaDrReal(i,j)  
           dRdSigmaLtd(i,j)=1./dSigmaDrLtd  
           SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)  
           SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)  
158    
159    #ifdef ALLOW_AUTODIFF_TAMC
160    CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
161    #endif
162    
163            DO j=1-Oly+1,sNy+Oly-1
164             DO i=1-Olx+1,sNx+Olx-1
165              IF ( dSigmaDrReal(i,j) .NE. 0. ) THEN
166               IF (dSigmaDrReal(i,j).GE.(-Small_Number))
167         &         dSigmaDrReal(i,j) = -Small_Number
168              ENDIF
169           ENDDO           ENDDO
170          ENDDO          ENDDO
171    
172        elseif (GM_taper_scheme.EQ.'gkw91') then  #ifdef ALLOW_AUTODIFF_TAMC
173    CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
174    CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
175    CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
176    #endif
177    
 C-      Gerdes, Koberle and Willebrand, Clim. Dyn. 1991  
178          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
179           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
180              IF ( dSigmaDrReal(i,j) .EQ. 0. ) THEN
181               IF ( SlopeX(i,j) .NE. 0. )
182         &      SlopeX(i,j) = SIGN(Small_taper,SlopeX(i,j))
183               IF ( SlopeY(i,j) .NE. 0. )
184         &      SlopeY(i,j) = SIGN(Small_taper,SlopeY(i,j))
185              ELSE
186               dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)
187               SlopeX(i,j) = -SlopeX(i,j)*dRdSigmaLtd(i,j)
188               SlopeY(i,j) = -SlopeY(i,j)*dRdSigmaLtd(i,j)
189              ENDIF
190             ENDDO
191            ENDDO
192    
193            gradSmod=SlopeX(i,j)*SlopeX(i,j)  #ifdef ALLOW_AUTODIFF_TAMC
194       &            +SlopeY(i,j)*SlopeY(i,j)  CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
195            if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)  CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
196    #endif
           dSigmaDrLtd=dSigmaDrReal(i,j)  
           if (dSigmaDrLtd.NE.0.) then  
            dRdSigmaLtd(i,j)=1./dSigmaDrLtd  
           else  
            dRdSigmaLtd(i,j)=0.  
           endif  
           if (gradSmod.LE.GM_maxSlope*abs(dSigmaDrReal(i,j))) then  
 C           If the slope is < Smax then calculate slopes properly  
             SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)  
             SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)  
           else  
 C           If the slope is > Smax then taper slopes  
             SlopeX(i,j)=-SlopeX(i,j)*dSigmaDrLtd  
      &                   *((GM_maxSlope/gradSmod)**2)  
             SlopeY(i,j)=-SlopeY(i,j)*dSigmaDrLtd  
      &                   *((GM_maxSlope/gradSmod)**2)  
           endif  
197    
198            DO j=1-Oly+1,sNy+Oly-1
199             DO i=1-Olx+1,sNx+Olx-1
200              SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
201         &                   +SlopeY(i,j)*SlopeY(i,j)
202              taperFct(i,j) = 1. _d 0
203           ENDDO           ENDDO
204          ENDDO          ENDDO
205    
206        elseif (GM_taper_scheme.EQ.'dm95') then  C- Compute the tapering function for the GM+Redi tensor :
207    
208  C-      Danabasoglu and McWilliams, J. Clim. 1995         IF (GM_taper_scheme.EQ.'linear') THEN
209    
210    C-      Simplest adiabatic tapering = Smax/Slope (linear)
211            maxSlopeSqr = GM_maxSlope*GM_maxSlope
212          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
213           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
214    
215            dSigmaDrLtd=dSigmaDrReal(i,j)            IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
216            if (dSigmaDrLtd.NE.0.) then             taperFct(i,j) = 1. _d 0
217             dRdSigmaLtd(i,j)=1./dSigmaDrLtd            ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr )  THEN
218            else             taperFct(i,j) = sqrt(maxSlopeSqr / SlopeSqr(i,j))
219             dRdSigmaLtd(i,j)=0.            ENDIF
           endif  
           sx   = -SlopeX(i,j)*dRdSigmaLtd(i,j)  
           sy   = -SlopeY(i,j)*dRdSigmaLtd(i,j)  
           Smod = sx*sx + sy*sy  
           if (Smod.NE.0.) Smod=sqrt(Smod)  
           f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))  
           SlopeX(i,j) = sx*f1  
           SlopeY(i,j) = sy*f1  
220    
221           ENDDO           ENDDO
222          ENDDO          ENDDO
223    
224        elseif (GM_taper_scheme.EQ.'ldd97') then         ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
225    
226  C-      Large, Danabasoglu and Doney, JPO 1997  C-      Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
227            maxSlopeSqr = GM_maxSlope*GM_maxSlope
228          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
229           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
230    
231            dSigmaDrLtd=dSigmaDrReal(i,j)            IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
232            if (dSigmaDrLtd.NE.0.) then             taperFct(i,j) = 1. _d 0
233             dRdSigmaLtd(i,j)=1./dSigmaDrLtd            ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr ) THEN
234            else             taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
235             dRdSigmaLtd(i,j)=0.            ENDIF
           endif  
           sx   = -SlopeX(i,j)*dRdSigmaLtd(i,j)  
           sy   = -SlopeY(i,j)*dRdSigmaLtd(i,j)  
           Smod = sx*sx + sy*sy  
           if (Smod.NE.0.) Smod=sqrt(Smod)  
           f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))  
           Cspd=2.  
           Lrho=100.e3  
           if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))  
           Lrho=min(Lrho,100.e3)  
           Lrho=max(Lrho,15.e3)  
           if (Smod.NE.0.) then  
             Rnondim=depthZ/(Lrho*Smod)  
           else  
             Rnondim=0.  
           endif  
           f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))  
           SlopeX(i,j) = sx*f1  
           SlopeY(i,j) = sy*f1  
236    
237           ENDDO           ENDDO
238          ENDDO          ENDDO
239    
240        elseif (GM_taper_scheme.EQ.' ') then         ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
241    
242  C-      No tapering/clipping selected so ...  C-      Danabasoglu and McWilliams, J. Clim. 1995
243            DO j=1-Oly+1,sNy+Oly-1
244             DO i=1-Olx+1,sNx+Olx-1
245    
246              IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
247               taperFct(i,j) = 1. _d 0
248              ELSE
249               Smod=sqrt(SlopeSqr(i,j))
250               taperFct(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
251              ENDIF
252             ENDDO
253            ENDDO
254    
255           ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
256    
257    C-      Large, Danabasoglu and Doney, JPO 1997
258          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
259           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
260    
261            dSigmaDrLtd=dSigmaDrReal(i,j)            IF (SlopeSqr(i,j) .EQ. 0.) THEN
262            if (dSigmaDrLtd.NE.0.) then             taperFct(i,j) = 1. _d 0
263              dRdSigmaLtd(i,j)=1./dSigmaDrLtd            ELSE
264            else             Smod=sqrt(SlopeSqr(i,j))
265              dRdSigmaLtd(i,j)=0.             f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
266            endif             Cspd=2.
267            SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)             Lrho=100. _d 3
268            SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)             if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))
269               Lrho=min(Lrho , 100. _d 3)
270               Lrho=max(Lrho , 15. _d 3)
271               Rnondim=depthZ/(Lrho*Smod)
272               f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))
273               taperFct(i,j)=f1*f2
274              ENDIF
275    
276           ENDDO           ENDDO
277          ENDDO          ENDDO
278    
279        endif         ELSEIF (GM_taper_scheme.NE.' ') THEN
280            STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'
281           ENDIF
282    
283          ENDIF
284    
285    
286  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22