/[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.11 by heimbach, Fri Feb 15 21:28:07 2002 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    
# Line 29  C     \================================= Line 27  C     \=================================
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"
# Line 52  CEndOfInterface Line 51  CEndOfInterface
51    
52  C     == Local variables ==  C     == Local variables ==
53        _RL Small_Number        _RL Small_Number
54          _RL Small_Taper
55        PARAMETER(Small_Number=1.D-12)        PARAMETER(Small_Number=1.D-12)
56          PARAMETER(Small_Taper=1.D+03)
57        _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58        _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59        _RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
# Line 63  C     == Local variables == Line 64  C     == Local variables ==
64        INTEGER i,j        INTEGER i,j
65    
66  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
67            act1 = bi - myBxLo(myThid)        act1 = bi - myBxLo(myThid)
68            max1 = myBxHi(myThid) - myBxLo(myThid) + 1        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
69            act2 = bj - myByLo(myThid)        act2 = bj - myByLo(myThid)
70            max2 = myByHi(myThid) - myByLo(myThid) + 1        max2 = myByHi(myThid) - myByLo(myThid) + 1
71            act3 = myThid - 1        act3 = myThid - 1
72            max3 = nTx*nTy        max3 = nTx*nTy
73            act4 = ikey_dynamics - 1        act4 = ikey_dynamics - 1
74            ikey = (act1 + 1) + act2*max1        ikey = (act1 + 1) + act2*max1
75       &                      + act3*max1*max2       &                  + act3*max1*max2
76       &                      + act4*max1*max2*max3       &                  + act4*max1*max2*max3
77  #endif /* ALLOW_AUTODIFF_TAMC */  #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.        IF (GM_taper_scheme.EQ.'orig' .OR.
87       &    GM_taper_scheme.EQ.'clipping') THEN       &    GM_taper_scheme.EQ.'clipping') THEN
88    
# Line 84  C       (this turns out to be the same a Line 92  C       (this turns out to be the same a
92  C-      Cox 1987 "Slope clipping"  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            gradSmod(i,j)=SlopeX(i,j)*SlopeX(i,j)            IF ( SlopeX(i,j)*SlopeX(i,j) + SlopeY(i,j)*SlopeY(i,j)
96       &            +SlopeY(i,j)*SlopeY(i,j)       &        .EQ. 0. ) THEN
97           ENDDO             gradSmod(i,j) = 0. _d 0
98          ENDDO            ELSE
99               gradSmod(i,j) = sqrt(SlopeX(i,j)*SlopeX(i,j)
100  #ifdef ALLOW_AUTODIFF_TAMC       &                         +SlopeY(i,j)*SlopeY(i,j))
101  CADJ STORE gradSmod(:,:)     = comlev1_bibj, key=ikey, byte=isbyte            ENDIF
 #endif  
   
         DO j=1-Oly+1,sNy+Oly-1  
          DO i=1-Olx+1,sNx+Olx-1  
           if (gradSmod(i,j) .NE. 0.) gradSmod(i,j)=sqrt(gradSmod(i,j))  
102           ENDDO           ENDDO
103          ENDDO          ENDDO
104    
# Line 106  CADJ STORE dSigmaDrReal(:,:) = comlev1_b Line 109  CADJ STORE dSigmaDrReal(:,:) = comlev1_b
109    
110          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
111           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
112            dSigmaDrLtd(i,j) = -(Small_Number+gradSmod(i,j)*GM_rMaxSlope)            IF (gradSmod(i,j) .NE. 0.) THEN
113            IF (dSigmaDrReal(i,j).GE.dSigmaDrLtd(i,j))             dSigmaDrLtd(i,j) = -gradSmod(i,j)*GM_rMaxSlope
114       &        dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)             IF (dSigmaDrReal(i,j) .GE. dSigmaDrLtd(i,j))
115         &         dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
116              ENDIF
117           ENDDO           ENDDO
118          ENDDO          ENDDO
119    
# Line 120  CADJ STORE dSigmaDrReal(:,:) = comlev1_b Line 125  CADJ STORE dSigmaDrReal(:,:) = comlev1_b
125    
126          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
127           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
128            dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)            IF (gradSmod(i,j) .EQ. 0.) THEN
129            SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)             SlopeX(i,j) = 0. _d 0
130            SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)             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    
# Line 136  CADJ STORE slopeY(:,:)       = comlev1_b Line 146  CADJ STORE slopeY(:,:)       = comlev1_b
146            SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)            SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)
147       &                 +SlopeY(i,j)*SlopeY(i,j)       &                 +SlopeY(i,j)*SlopeY(i,j)
148            taperFct(i,j)=1. _d 0            taperFct(i,j)=1. _d 0
   
149           ENDDO           ENDDO
150          ENDDO          ENDDO
151    
152        ELSE        ELSE
153    
154    C----------------------------------------------------------------------
155    
156  C- Compute the slope, no clipping, but avoid reverse slope in negatively  C- Compute the slope, no clipping, but avoid reverse slope in negatively
157  C                                  stratified (Sigma_Z > 0) region :  C                                  stratified (Sigma_Z > 0) region :
158    
# Line 151  CADJ STORE dSigmaDrReal(:,:) = comlev1_b Line 162  CADJ STORE dSigmaDrReal(:,:) = comlev1_b
162    
163          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
164           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
165            dSigmaDrLtd(i,j) = -Small_Number            IF ( dSigmaDrReal(i,j) .NE. 0. ) THEN
166            IF (dSigmaDrReal(i,j).GE.dSigmaDrLtd(i,j))             IF (dSigmaDrReal(i,j).GE.(-Small_Number))
167       &        dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)       &         dSigmaDrReal(i,j) = -Small_Number
168              ENDIF
169           ENDDO           ENDDO
170          ENDDO          ENDDO
171    
# Line 165  CADJ STORE dSigmaDrReal(:,:) = comlev1_b Line 177  CADJ STORE dSigmaDrReal(:,:) = comlev1_b
177    
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            dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)            IF ( dSigmaDrReal(i,j) .EQ. 0. ) THEN
181            SlopeX(i,j) = -SlopeX(i,j)*dRdSigmaLtd(i,j)             IF ( SlopeX(i,j) .NE. 0. )
182            SlopeY(i,j) = -SlopeY(i,j)*dRdSigmaLtd(i,j)       &      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           ENDDO
191          ENDDO          ENDDO
192    
# Line 178  CADJ STORE slopeY(:,:)       = comlev1_b Line 197  CADJ STORE slopeY(:,:)       = comlev1_b
197    
198          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
199           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
200            SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)            SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
201       &                 +SlopeY(i,j)*SlopeY(i,j)       &                   +SlopeY(i,j)*SlopeY(i,j)
202            taperFct(i,j)=1. _d 0            taperFct(i,j) = 1. _d 0
203           ENDDO           ENDDO
204          ENDDO          ENDDO
205    
# Line 193  C-      Simplest adiabatic tapering = Sm Line 212  C-      Simplest adiabatic tapering = Sm
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            IF (SlopeSqr(i,j).GT.maxSlopeSqr)            IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
216       &      taperFct(i,j)=GM_maxSlope/sqrt(SlopeSqr(i,j))             taperFct(i,j) = 1. _d 0
217              ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr )  THEN
218               taperFct(i,j) = sqrt(maxSlopeSqr / SlopeSqr(i,j))
219              ENDIF
220    
221           ENDDO           ENDDO
222          ENDDO          ENDDO
# Line 206  C-      Gerdes, Koberle and Willebrand, Line 228  C-      Gerdes, Koberle and Willebrand,
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            IF (SlopeSqr(i,j).GT.maxSlopeSqr)            IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
232       &      taperFct(i,j)=maxSlopeSqr/SlopeSqr(i,j)             taperFct(i,j) = 1. _d 0
233              ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr ) THEN
234               taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
235              ENDIF
236    
237           ENDDO           ENDDO
238          ENDDO          ENDDO
# Line 218  C-      Danabasoglu and McWilliams, J. C Line 243  C-      Danabasoglu and McWilliams, J. C
243          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
244           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
245    
246            Smod = SlopeSqr(i,j)            IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
247            if (Smod.NE.0.) Smod=sqrt(Smod)             taperFct(i,j) = 1. _d 0
248            taperFct(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))            ELSE
249               Smod=sqrt(SlopeSqr(i,j))
250               taperFct(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
251              ENDIF
252           ENDDO           ENDDO
253          ENDDO          ENDDO
254    
# Line 231  C-      Large, Danabasoglu and Doney, JP Line 258  C-      Large, Danabasoglu and Doney, JP
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            Smod = SlopeSqr(i,j)            IF (SlopeSqr(i,j) .EQ. 0.) THEN
262            if (Smod.NE.0.) Smod=sqrt(Smod)             taperFct(i,j) = 1. _d 0
263            f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))            ELSE
264            Cspd=2.             Smod=sqrt(SlopeSqr(i,j))
265            Lrho=100.e3             f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
266            if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))             Cspd=2.
267            Lrho=min(Lrho , 100. _d 3)             Lrho=100. _d 3
268            Lrho=max(Lrho , 15. _d 3)             if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))
269            if (Smod.NE.0.) then             Lrho=min(Lrho , 100. _d 3)
270              Rnondim=depthZ/(Lrho*Smod)             Lrho=max(Lrho , 15. _d 3)
271            else             Rnondim=depthZ/(Lrho*Smod)
272              Rnondim=0.             f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))
273            endif             taperFct(i,j)=f1*f2
274            f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))            ENDIF
           taperFct(i,j)=f1*f2  
275    
276           ENDDO           ENDDO
277          ENDDO          ENDDO
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22