/[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.5 by jmc, Sun Jan 12 21:27:20 2003 UTC revision 1.6 by heimbach, Tue Jan 21 19:34:13 2003 UTC
# Line 56  C     == Local variables == Line 56  C     == Local variables ==
56        _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57        _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
58        _RL drdsigmaltd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL drdsigmaltd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59          _RL  SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60        _RL f1,Smod,f2,Rnondim,Cspd,Lrho        _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
# Line 70  C     == Local variables == Line 72  C     == Local variables ==
72        act3 = myThid - 1        act3 = myThid - 1
73        max3 = nTx*nTy        max3 = nTx*nTy
74        act4 = ikey_dynamics - 1        act4 = ikey_dynamics - 1
75        ikey = (act1 + 1) + act2*max1        igmkey = (act1 + 1) + act2*max1
76       &                  + act3*max1*max2       &                    + act3*max1*max2
77       &                  + act4*max1*max2*max3       &                    + act4*max1*max2*max3
78        kkey = (ikey-1)*Nr + k        kkey = (igmkey-1)*Nr + k
79  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
80    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE slopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif  
   
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    
# Line 159  CADJ STORE dSigmaDrS(:,:)    = comlev1_b Line 162  CADJ STORE dSigmaDrS(:,:)    = comlev1_b
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  #ifdef ALLOW_AUTODIFF_TAMC
176  CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE slopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
177  CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDrW(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
178  #endif  #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
# Line 182  CADJ STORE dsigmadrW(:,:)    = comlev1_b Line 195  CADJ STORE dsigmadrW(:,:)    = comlev1_b
195           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
196            drdsigmaltd(i,j) = 1./dSigmaDrW(i,j)            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            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           ENDDO
211          ENDDO          ENDDO
212    
# Line 207  CADJ STORE dsigmadrS(:,:)    = comlev1_b Line 231  CADJ STORE dsigmadrS(:,:)    = comlev1_b
231            taperY(i,j)=1. _d 0            taperY(i,j)=1. _d 0
232           ENDDO           ENDDO
233          ENDDO          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
244            ENDDO
245    
246  C- Compute the tapering function for the GM+Redi tensor :  C- Compute the tapering function for the GM+Redi tensor :
247    
# Line 221  C-      Simplest adiabatic tapering = Sm Line 256  C-      Simplest adiabatic tapering = Sm
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 232  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 271  C-      Large, Danabasoglu and Doney, JP Line 315  C-      Large, Danabasoglu and Doney, JP
315            Lrho=max(Lrho , 15. _d 3)            Lrho=max(Lrho , 15. _d 3)
316    
317            Smod = abs(SlopeX(i,j))            Smod = abs(SlopeX(i,j))
318              if ( Smod .LT. tmpvar ) then
319            f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))            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)
# Line 279  C-      Large, Danabasoglu and Doney, JP Line 324  C-      Large, Danabasoglu and Doney, JP
324            endif            endif
325            f2=op5*( 1. _d 0 + sin( fpi*(Rnondim-op5)))            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              if ( Smod .LT. tmpvar ) then
331            f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))            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)
# Line 289  C-      Large, Danabasoglu and Doney, JP Line 336  C-      Large, Danabasoglu and Doney, JP
336            endif            endif
337            f2=op5*( 1. _d 0 + sin( fpi*(Rnondim-op5)))            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 297  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    
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22