pkg/gmredi/gmredi_slope_limit.F

C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_slope_limit.F,v 1.9 2001/12/16 18:54:49 jmc Exp $
C $Name:  $

#include "GMREDI_OPTIONS.h"

CStartOfInterface
      SUBROUTINE GMREDI_SLOPE_LIMIT(
     I             dSigmaDrReal,
     I             depthZ,
     U             SlopeX, SlopeY,
     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 "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 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,myThid
CEndOfInterface

#ifdef ALLOW_GMREDI

C     == Local variables ==
      _RL Small_Number
      PARAMETER(Small_Number=1.D-12)
      _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 */

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

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
          gradSmod(i,j)=SlopeX(i,j)*SlopeX(i,j)
     &            +SlopeY(i,j)*SlopeY(i,j)
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE gradSmod(:,:)     = 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.) gradSmod(i,j)=sqrt(gradSmod(i,j))
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE gradSmod(:,:)     = comlev1_bibj, key=ikey, byte=isbyte
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
          dSigmaDrLtd(i,j) = -(Small_Number+gradSmod(i,j)*GM_rMaxSlope)
          IF (dSigmaDrReal(i,j).GE.dSigmaDrLtd(i,j))
     &        dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
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
          dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)
          SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)
          SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
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

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

#ifdef ALLOW_AUTODIFF_TAMC
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
          dSigmaDrLtd(i,j) = -Small_Number
          IF (dSigmaDrReal(i,j).GE.dSigmaDrLtd(i,j))
     &        dSigmaDrReal(i,j) = dSigmaDrLtd(i,j)
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE slopeY(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
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
          dRdSigmaLtd(i,j) = 1./dSigmaDrReal(i,j)
          SlopeX(i,j) = -SlopeX(i,j)*dRdSigmaLtd(i,j)
          SlopeY(i,j) = -SlopeY(i,j)*dRdSigmaLtd(i,j)
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE slopeX(:,:)       = comlev1_bibj, key=ikey, byte=isbyte
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).GT.maxSlopeSqr) 
     &      taperFct(i,j)=GM_maxSlope/sqrt(SlopeSqr(i,j))

         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).GT.maxSlopeSqr)
     &      taperFct(i,j)=maxSlopeSqr/SlopeSqr(i,j)

         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

          Smod = SlopeSqr(i,j)
          if (Smod.NE.0.) Smod=sqrt(Smod)
          taperFct(i,j)=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))

         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

          Smod = SlopeSqr(i,j)
          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)))
          taperFct(i,j)=f1*f2

         ENDDO
        ENDDO

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

      ENDIF


#endif /* ALLOW_GMREDI */

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