pkg/gmredi/gmredi_slope_limit.F


#include "GMREDI_OPTIONS.h"

CStartOfInterface
      SUBROUTINE GMREDI_SLOPE_LIMIT(
     I             dSigmaDrReal,
     I             depthZ,K,
     U             SlopeX, SlopeY,
     U             dSigmaDx, dSigmaDy,
     O             SlopeSqr, taperFct,
     I             bi,bj, myThid )
C     /==========================================================\
C     | SUBROUTINE GMREDI_SLOPE_LIMIT                            |
C     | o Calculate slopes for use in GM/Redi tensor             |
C     |==========================================================|
C     | On entry:                                                |
C     |            dSigmaDrReal conatins the d/dz Sigma          |
C     |            SlopeX/Y     contains X/Y gradients of sigma  |
C     |            depthZ       conatins the height (m) of level |
C     | On exit:                                                 |
C     |            dSigmaDrReal conatins the effective dSig/dz   |
C     |            SlopeX/Y     contains X/Y slopes              |
C     |            SlopeSqr     contains Sx^2+Sy^2               |
C     |            taperFct     contains tapering funct. value ; |
C     |                         = 1 when using no tapering       |
C     \==========================================================/
      IMPLICIT NONE

C     == Global variables ==
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "GMREDI.h"
#include "PARAMS.h"
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */

C     == Routine arguments ==
C
      _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL depthZ
      INTEGER bi,bj,K,myThid
CEndOfInterface

#ifdef ALLOW_GMREDI

C     == Local variables ==
      _RL Small_Number
      _RL Small_Taper
      PARAMETER(Small_Number=1.D-12)
      PARAMETER(Small_Taper=1.D+03)
      _RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL f1,Smod,f2,Rnondim,Cspd,Lrho
      _RL maxSlopeSqr
      _RL fpi
      PARAMETER(fpi=3.141592653589793047592d0)
      INTEGER i,j

#ifdef ALLOW_AUTODIFF_TAMC
      act1 = bi - myBxLo(myThid)
      max1 = myBxHi(myThid) - myBxLo(myThid) + 1
      act2 = bj - myByLo(myThid)
      max2 = myByHi(myThid) - myByLo(myThid) + 1
      act3 = myThid - 1
      max3 = nTx*nTy
      act4 = ikey_dynamics - 1
      ikey = (act1 + 1) + act2*max1
     &                  + act3*max1*max2
     &                  + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

      DO j=1-Oly+1,sNy+Oly-1
       DO i=1-Olx+1,sNx+Olx-1
        gradSmod(i,j)    = 0. _d 0
        dSigmaDrLtd(i,j) = 0. _d 0
       ENDDO
      ENDDO

      IF (GM_taper_scheme.EQ.'orig' .OR.
     &    GM_taper_scheme.EQ.'clipping') THEN

#if (defined (GM_TAPER_ORIG_CLIPPING))

C-      Original implementation in mitgcmuv
C       (this turns out to be the same as Cox slope clipping)

C-      Cox 1987 "Slope clipping"
        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1
          IF ( DSigmaDX(i,j)*DSigmaDX(i,j) + 
     &        DSigmaDY(i,j)*DSigmaDY(i,j)
     &        .EQ. 0. ) THEN
           gradSmod(i,j) = 0. _d 0
          ELSE
           gradSmod(i,j) = sqrt(DSigmaDX(i,j)*DSigmaDX(i,j)
     &                         +DSigmaDY(i,j)*DSigmaDY(i,j))
          ENDIF
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE gradSmod(:,:)     = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif

        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1
          IF (gradSmod(i,j) .NE. 0.) THEN
           dSigmaDrLtd(i,j) = -gradSmod(i,j)*GM_rMaxSlope
           IF ( dSigmaDrReal(i,j) .GE. 
     &          GM_adjointRescale*dSigmaDrLtd(i,j) )
     &          dSigmaDrReal(i,j) = 
     &          dSigmaDrLtd(i,j)*GM_adjointRescale
ctest           dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
          ENDIF
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif

        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1
          IF (gradSmod(i,j) .EQ. 0.) THEN
           SlopeX(i,j) = 0. _d 0
           SlopeY(i,j) = 0. _d 0
          ELSE
           dRdSigmaLtd(i,j) = 1./( dSigmaDrReal(i,j) )
           SlopeX(i,j)=-DSigmaDX(i,j)*dRdSigmaLtd(i,j)
           SlopeY(i,j)=-DSigmaDY(i,j)*dRdSigmaLtd(i,j)
          ENDIF
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
#endif

        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1
          SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)
     &                 +SlopeY(i,j)*SlopeY(i,j)
          taperFct(i,j)=1. _d 0
         ENDDO
        ENDDO

#else

        STOP 'Need to set runtime flag GM_taper_scheme correctly'

#endif /* GM_TAPER_ORIG_CLIPPING */

      ELSE IF (GM_taper_scheme.EQ.'ac02') THEN

#if (defined (GM_TAPER_AC02))

        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1
          dRdSigmaLtd(i,j)= 
     &                      dSigmaDx(i,j)*dSigmaDx(i,j)
     &                      + dSigmaDy(i,j)*dSigmaDy(i,j)
     &                      + dSigmaDrReal(i,j)**2
ctest
          IF (dRdSigmaLtd(i,j).NE.0.) then
             dRdSigmaLtd(i,j)=1. _d 0
     &            / ( dRdSigmaLtd(i,j) )
             SlopeSqr(i,j)=(dSigmaDx(i,j)*dSigmaDx(i,j)
     &            +dSigmaDy(i,j)*dSigmaDy(i,j))*dRdSigmaLtd(i,j)
             SlopeX(i,j)=-dSigmaDx(i,j)
     &            *dRdSigmaLtd(i,j)*dSigmaDrReal(i,j)
             SlopeY(i,j)=-dSigmaDy(i,j)
     &            *dRdSigmaLtd(i,j)*dSigmaDrReal(i,j)
cph             T11(i,j)=(dSigmaDrReal(i,j)**2)*dRdSigmaLtd(i,j)
             taperFct(i,j) = 1. _d 0
cph          ELSE
cph             SlopeSqr(i,j) = 0. _d 0
cph             SlopeX(i,j)   = 0. _d 0
cph             SlopeY(i,j)   = 0. _d 0
cph             T11(i,j)      = 0. _d 0
cph             taperFct(i,j) = 0. _d 0
          ENDIF
cph          IF ( SlopeSqr(i,j) .GT. GM_maxSlope**2 ) THEN
cph           taperFct(i,j) = GM_maxSlope**2/SlopeSqr(i,j)
cph          ENDIF
         ENDDO
        ENDDO

#else
       
        STOP 'Need to set runtime flag GM_taper_scheme correctly'

#endif /* GM_TAPER_AC02 */

      ELSE

#if (defined (GM_TAPER_REST))

C----------------------------------------------------------------------

C- Compute the slope, no clipping, but avoid reverse slope in negatively
C                                  stratified (Sigma_Z > 0) region :

#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif

        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1
          IF ( dSigmaDrReal(i,j) .NE. 0. ) THEN
           IF (dSigmaDrReal(i,j).GE.(-Small_Number))
     &         dSigmaDrReal(i,j) = -Small_Number
          ENDIF
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE dSigmaDx(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE dSigmaDy(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif

        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1
          IF ( dSigmaDrReal(i,j) .EQ. 0. ) THEN
           IF ( dSigmaDx(i,j) .NE. 0. )
     &      SlopeX(i,j) = SIGN(Small_taper,dSigmaDx(i,j))
           IF ( dSigmaDy(i,j) .NE. 0. )
     &      SlopeY(i,j) = SIGN(Small_taper,dSigmaDy(i,j))
          ELSE
           dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)
           SlopeX(i,j) = -dSigmaDx(i,j)*dRdSigmaLtd(i,j)
           SlopeY(i,j) = -dSigmaDy(i,j)*dRdSigmaLtd(i,j)
          ENDIF
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
#endif

        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1
          SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
     &                   +SlopeY(i,j)*SlopeY(i,j)
          taperFct(i,j) = 1. _d 0
         ENDDO
        ENDDO

C- Compute the tapering function for the GM+Redi tensor :

       IF (GM_taper_scheme.EQ.'linear') THEN

C-      Simplest adiabatic tapering = Smax/Slope (linear)
        maxSlopeSqr = GM_maxSlope*GM_maxSlope
        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1

          IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
           taperFct(i,j) = 1. _d 0
          ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr )  THEN
           taperFct(i,j) = sqrt(maxSlopeSqr / SlopeSqr(i,j))
          ENDIF

         ENDDO
        ENDDO

       ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN

C-      Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
        maxSlopeSqr = GM_maxSlope*GM_maxSlope
        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1

          IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
           taperFct(i,j) = 1. _d 0
          ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr ) THEN
           taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
          ENDIF

         ENDDO
        ENDDO

       ELSEIF (GM_taper_scheme.EQ.'dm95') THEN

C-      Danabasoglu and McWilliams, J. Clim. 1995
        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1

          IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
           taperFct(i,j) = 1. _d 0
          ELSE
           Smod=sqrt(SlopeSqr(i,j))
           taperFct(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
          ENDIF
         ENDDO
        ENDDO

       ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN

C-      Large, Danabasoglu and Doney, JPO 1997
        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1

          IF (SlopeSqr(i,j) .EQ. 0.) THEN
           taperFct(i,j) = 1. _d 0
          ELSE
           Smod=sqrt(SlopeSqr(i,j))
           f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
           Cspd=2.
           Lrho=100. _d 3
           if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))
           Lrho=min(Lrho , 100. _d 3)
           Lrho=max(Lrho , 15. _d 3)
           Rnondim=depthZ/(Lrho*Smod)
           f2=0.5*(1.+sin( fpi*(Rnondim-0.5)))
           taperFct(i,j)=f1*f2
          ENDIF

         ENDDO
        ENDDO

       ELSEIF (GM_taper_scheme.NE.' ') THEN
        STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'
       ENDIF

#endif

      ENDIF


#endif /* ALLOW_GMREDI */

      RETURN
      END
1	adcroft	1.1
2			#include "GMREDI_OPTIONS.h"
3
4			CStartOfInterface
5			SUBROUTINE GMREDI_SLOPE_LIMIT(
6	jmc	1.8	I dSigmaDrReal,
7	heimbach	1.13	I depthZ,K,
8	adcroft	1.1	U SlopeX, SlopeY,
9	heimbach	1.13	U dSigmaDx, dSigmaDy,
10	jmc	1.8	O SlopeSqr, taperFct,
11	adcroft	1.1	I bi,bj, myThid )
12			C /==========================================================\
13			C \| SUBROUTINE GMREDI_SLOPE_LIMIT \|
14			C \| o Calculate slopes for use in GM/Redi tensor \|
15			C \|==========================================================\|
16			C \| On entry: \|
17	jmc	1.8	C \| dSigmaDrReal conatins the d/dz Sigma \|
18	adcroft	1.1	C \| SlopeX/Y contains X/Y gradients of sigma \|
19			C \| depthZ conatins the height (m) of level \|
20			C \| On exit: \|
21	jmc	1.8	C \| dSigmaDrReal conatins the effective dSig/dz \|
22	adcroft	1.1	C \| SlopeX/Y contains X/Y slopes \|
23	jmc	1.8	C \| SlopeSqr contains Sx^2+Sy^2 \|
24			C \| taperFct contains tapering funct. value ; \|
25			C \| = 1 when using no tapering \|
26	adcroft	1.1	C \==========================================================/
27			IMPLICIT NONE
28
29			C == Global variables ==
30			#include "SIZE.h"
31	heimbach	1.12	#include "GRID.h"
32	adcroft	1.1	#include "EEPARAMS.h"
33			#include "GMREDI.h"
34			#include "PARAMS.h"
35	heimbach	1.10	#ifdef ALLOW_AUTODIFF_TAMC
36			#include "tamc.h"
37			#include "tamc_keys.h"
38			#endif /* ALLOW_AUTODIFF_TAMC */
39	adcroft	1.1
40			C == Routine arguments ==
41			C
42			_RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43			_RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44	heimbach	1.13	_RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45			_RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46	adcroft	1.1	_RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47	jmc	1.8	_RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48			_RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49	adcroft	1.1	_RL depthZ
50	heimbach	1.13	INTEGER bi,bj,K,myThid
51	adcroft	1.1	CEndOfInterface
52
53			#ifdef ALLOW_GMREDI
54
55			C == Local variables ==
56			_RL Small_Number
57	heimbach	1.12	_RL Small_Taper
58	heimbach	1.7	PARAMETER(Small_Number=1.D-12)
59	heimbach	1.12	PARAMETER(Small_Taper=1.D+03)
60	heimbach	1.10	_RL gradSmod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
61			_RL dSigmaDrLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
62			_RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
63			_RL f1,Smod,f2,Rnondim,Cspd,Lrho
64	jmc	1.8	_RL maxSlopeSqr
65	adcroft	1.1	_RL fpi
66	heimbach	1.6	PARAMETER(fpi=3.141592653589793047592d0)
67	adcroft	1.1	INTEGER i,j
68
69	heimbach	1.10	#ifdef ALLOW_AUTODIFF_TAMC
70	heimbach	1.12	act1 = bi - myBxLo(myThid)
71			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
72			act2 = bj - myByLo(myThid)
73			max2 = myByHi(myThid) - myByLo(myThid) + 1
74			act3 = myThid - 1
75			max3 = nTx*nTy
76			act4 = ikey_dynamics - 1
77			ikey = (act1 + 1) + act2*max1
78			& + act3max1max2
79			& + act4max1max2*max3
80	heimbach	1.10	#endif /* ALLOW_AUTODIFF_TAMC */
81
82	heimbach	1.12	DO j=1-Oly+1,sNy+Oly-1
83			DO i=1-Olx+1,sNx+Olx-1
84			gradSmod(i,j) = 0. _d 0
85			dSigmaDrLtd(i,j) = 0. _d 0
86			ENDDO
87			ENDDO
88
89	jmc	1.8	IF (GM_taper_scheme.EQ.'orig' .OR.
90			& GM_taper_scheme.EQ.'clipping') THEN
91	adcroft	1.1
92	heimbach	1.13	#if (defined (GM_TAPER_ORIG_CLIPPING))
93
94	adcroft	1.1	C- Original implementation in mitgcmuv
95			C (this turns out to be the same as Cox slope clipping)
96	jmc	1.8
97			C- Cox 1987 "Slope clipping"
98	adcroft	1.1	DO j=1-Oly+1,sNy+Oly-1
99			DO i=1-Olx+1,sNx+Olx-1
100	heimbach	1.13	IF ( DSigmaDX(i,j)*DSigmaDX(i,j) +
101			& DSigmaDY(i,j)*DSigmaDY(i,j)
102	heimbach	1.12	& .EQ. 0. ) THEN
103			gradSmod(i,j) = 0. _d 0
104			ELSE
105	heimbach	1.13	gradSmod(i,j) = sqrt(DSigmaDX(i,j)*DSigmaDX(i,j)
106			& +DSigmaDY(i,j)*DSigmaDY(i,j))
107	heimbach	1.12	ENDIF
108	heimbach	1.10	ENDDO
109			ENDDO
110	adcroft	1.1
111	heimbach	1.10	#ifdef ALLOW_AUTODIFF_TAMC
112	heimbach	1.13	cnostore CADJ STORE gradSmod(:,:) = comlev1_bibj, key=ikey, byte=isbyte
113			cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
114	heimbach	1.10	#endif
115
116			DO j=1-Oly+1,sNy+Oly-1
117			DO i=1-Olx+1,sNx+Olx-1
118	heimbach	1.12	IF (gradSmod(i,j) .NE. 0.) THEN
119			dSigmaDrLtd(i,j) = -gradSmod(i,j)*GM_rMaxSlope
120	heimbach	1.13	IF ( dSigmaDrReal(i,j) .GE.
121			& GM_adjointRescale*dSigmaDrLtd(i,j) )
122			& dSigmaDrReal(i,j) =
123			& dSigmaDrLtd(i,j)*GM_adjointRescale
124			ctest dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
125	heimbach	1.12	ENDIF
126	heimbach	1.10	ENDDO
127			ENDDO
128
129			#ifdef ALLOW_AUTODIFF_TAMC
130	heimbach	1.13	cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
131			cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
132			cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
133	heimbach	1.10	#endif
134
135			DO j=1-Oly+1,sNy+Oly-1
136			DO i=1-Olx+1,sNx+Olx-1
137	heimbach	1.12	IF (gradSmod(i,j) .EQ. 0.) THEN
138			SlopeX(i,j) = 0. _d 0
139			SlopeY(i,j) = 0. _d 0
140			ELSE
141	heimbach	1.13	dRdSigmaLtd(i,j) = 1./( dSigmaDrReal(i,j) )
142			SlopeX(i,j)=-DSigmaDX(i,j)*dRdSigmaLtd(i,j)
143			SlopeY(i,j)=-DSigmaDY(i,j)*dRdSigmaLtd(i,j)
144	heimbach	1.12	ENDIF
145	heimbach	1.10	ENDDO
146			ENDDO
147	adcroft	1.1
148	heimbach	1.10	#ifdef ALLOW_AUTODIFF_TAMC
149	heimbach	1.13	cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
150			cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
151	heimbach	1.10	#endif
152	jmc	1.8
153	heimbach	1.10	DO j=1-Oly+1,sNy+Oly-1
154			DO i=1-Olx+1,sNx+Olx-1
155	jmc	1.8	SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)
156			& +SlopeY(i,j)*SlopeY(i,j)
157			taperFct(i,j)=1. _d 0
158	adcroft	1.1	ENDDO
159			ENDDO
160
161	heimbach	1.13	#else
162
163			STOP 'Need to set runtime flag GM_taper_scheme correctly'
164
165			#endif /* GM_TAPER_ORIG_CLIPPING */
166
167			ELSE IF (GM_taper_scheme.EQ.'ac02') THEN
168
169			#if (defined (GM_TAPER_AC02))
170
171			DO j=1-Oly+1,sNy+Oly-1
172			DO i=1-Olx+1,sNx+Olx-1
173			dRdSigmaLtd(i,j)=
174			& dSigmaDx(i,j)*dSigmaDx(i,j)
175			& + dSigmaDy(i,j)*dSigmaDy(i,j)
176			& + dSigmaDrReal(i,j)**2
177			ctest
178			IF (dRdSigmaLtd(i,j).NE.0.) then
179			dRdSigmaLtd(i,j)=1. _d 0
180			& / ( dRdSigmaLtd(i,j) )
181			SlopeSqr(i,j)=(dSigmaDx(i,j)*dSigmaDx(i,j)
182			& +dSigmaDy(i,j)dSigmaDy(i,j))dRdSigmaLtd(i,j)
183			SlopeX(i,j)=-dSigmaDx(i,j)
184			& dRdSigmaLtd(i,j)dSigmaDrReal(i,j)
185			SlopeY(i,j)=-dSigmaDy(i,j)
186			& dRdSigmaLtd(i,j)dSigmaDrReal(i,j)
187			cph T11(i,j)=(dSigmaDrReal(i,j)*2)dRdSigmaLtd(i,j)
188			taperFct(i,j) = 1. _d 0
189			cph ELSE
190			cph SlopeSqr(i,j) = 0. _d 0
191			cph SlopeX(i,j) = 0. _d 0
192			cph SlopeY(i,j) = 0. _d 0
193			cph T11(i,j) = 0. _d 0
194			cph taperFct(i,j) = 0. _d 0
195			ENDIF
196			cph IF ( SlopeSqr(i,j) .GT. GM_maxSlope**2 ) THEN
197			cph taperFct(i,j) = GM_maxSlope**2/SlopeSqr(i,j)
198			cph ENDIF
199			ENDDO
200			ENDDO
201
202			#else
203
204			STOP 'Need to set runtime flag GM_taper_scheme correctly'
205
206			#endif /* GM_TAPER_AC02 */
207
208	jmc	1.8	ELSE
209	adcroft	1.1
210	heimbach	1.13	#if (defined (GM_TAPER_REST))
211
212	heimbach	1.12	C----------------------------------------------------------------------
213
214	heimbach	1.10	C- Compute the slope, no clipping, but avoid reverse slope in negatively
215	jmc	1.8	C stratified (Sigma_Z > 0) region :
216	heimbach	1.10
217			#ifdef ALLOW_AUTODIFF_TAMC
218	heimbach	1.13	cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
219	heimbach	1.10	#endif
220
221			DO j=1-Oly+1,sNy+Oly-1
222			DO i=1-Olx+1,sNx+Olx-1
223	heimbach	1.12	IF ( dSigmaDrReal(i,j) .NE. 0. ) THEN
224			IF (dSigmaDrReal(i,j).GE.(-Small_Number))
225			& dSigmaDrReal(i,j) = -Small_Number
226			ENDIF
227	heimbach	1.10	ENDDO
228			ENDDO
229
230			#ifdef ALLOW_AUTODIFF_TAMC
231	heimbach	1.13	cnostore CADJ STORE dSigmaDx(:,:) = comlev1_bibj, key=ikey, byte=isbyte
232			cnostore CADJ STORE dSigmaDy(:,:) = comlev1_bibj, key=ikey, byte=isbyte
233			cnostore CADJ STORE dSigmaDrReal(:,:) = comlev1_bibj, key=ikey, byte=isbyte
234	heimbach	1.10	#endif
235
236	adcroft	1.1	DO j=1-Oly+1,sNy+Oly-1
237			DO i=1-Olx+1,sNx+Olx-1
238	heimbach	1.12	IF ( dSigmaDrReal(i,j) .EQ. 0. ) THEN
239	heimbach	1.13	IF ( dSigmaDx(i,j) .NE. 0. )
240			& SlopeX(i,j) = SIGN(Small_taper,dSigmaDx(i,j))
241			IF ( dSigmaDy(i,j) .NE. 0. )
242			& SlopeY(i,j) = SIGN(Small_taper,dSigmaDy(i,j))
243	heimbach	1.12	ELSE
244			dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)
245	heimbach	1.13	SlopeX(i,j) = -dSigmaDx(i,j)*dRdSigmaLtd(i,j)
246			SlopeY(i,j) = -dSigmaDy(i,j)*dRdSigmaLtd(i,j)
247	heimbach	1.12	ENDIF
248	heimbach	1.10	ENDDO
249			ENDDO
250	adcroft	1.1
251	heimbach	1.10	#ifdef ALLOW_AUTODIFF_TAMC
252	heimbach	1.13	cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
253			cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
254	heimbach	1.10	#endif
255	jmc	1.8
256	heimbach	1.10	DO j=1-Oly+1,sNy+Oly-1
257			DO i=1-Olx+1,sNx+Olx-1
258	heimbach	1.12	SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
259			& +SlopeY(i,j)*SlopeY(i,j)
260			taperFct(i,j) = 1. _d 0
261	jmc	1.8	ENDDO
262			ENDDO
263
264			C- Compute the tapering function for the GM+Redi tensor :
265
266			IF (GM_taper_scheme.EQ.'linear') THEN
267
268			C- Simplest adiabatic tapering = Smax/Slope (linear)
269			maxSlopeSqr = GM_maxSlope*GM_maxSlope
270			DO j=1-Oly+1,sNy+Oly-1
271			DO i=1-Olx+1,sNx+Olx-1
272	adcroft	1.1
273	heimbach	1.12	IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
274			taperFct(i,j) = 1. _d 0
275			ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr ) THEN
276			taperFct(i,j) = sqrt(maxSlopeSqr / SlopeSqr(i,j))
277			ENDIF
278	adcroft	1.1
279			ENDDO
280			ENDDO
281
282	jmc	1.8	ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
283	adcroft	1.1
284			C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
285	jmc	1.8	maxSlopeSqr = GM_maxSlope*GM_maxSlope
286	adcroft	1.1	DO j=1-Oly+1,sNy+Oly-1
287			DO i=1-Olx+1,sNx+Olx-1
288
289	heimbach	1.12	IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
290			taperFct(i,j) = 1. _d 0
291			ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr ) THEN
292			taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
293			ENDIF
294	adcroft	1.1
295			ENDDO
296			ENDDO
297
298	jmc	1.8	ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
299	adcroft	1.1
300			C- Danabasoglu and McWilliams, J. Clim. 1995
301			DO j=1-Oly+1,sNy+Oly-1
302			DO i=1-Olx+1,sNx+Olx-1
303
304	heimbach	1.12	IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
305			taperFct(i,j) = 1. _d 0
306			ELSE
307			Smod=sqrt(SlopeSqr(i,j))
308			taperFct(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
309			ENDIF
310	adcroft	1.1	ENDDO
311			ENDDO
312
313	jmc	1.8	ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
314	adcroft	1.1
315			C- Large, Danabasoglu and Doney, JPO 1997
316			DO j=1-Oly+1,sNy+Oly-1
317			DO i=1-Olx+1,sNx+Olx-1
318
319	heimbach	1.12	IF (SlopeSqr(i,j) .EQ. 0.) THEN
320			taperFct(i,j) = 1. _d 0
321			ELSE
322			Smod=sqrt(SlopeSqr(i,j))
323			f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
324			Cspd=2.
325			Lrho=100. _d 3
326			if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))
327			Lrho=min(Lrho , 100. _d 3)
328			Lrho=max(Lrho , 15. _d 3)
329			Rnondim=depthZ/(Lrho*Smod)
330			f2=0.5(1.+sin( fpi(Rnondim-0.5)))
331			taperFct(i,j)=f1*f2
332			ENDIF
333	adcroft	1.1
334			ENDDO
335			ENDDO
336
337	jmc	1.9	ELSEIF (GM_taper_scheme.NE.' ') THEN
338			STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'
339	jmc	1.8	ENDIF
340	heimbach	1.13
341			#endif
342	adcroft	1.1
343	jmc	1.8	ENDIF
344	adcroft	1.1
345
346			#endif /* ALLOW_GMREDI */
347
348			RETURN
349			END