pkg/gmredi/gmredi_slope_limit.F

C $Header: /u/gcmpack/models/MITgcmUV/pkg/gmredi/gmredi_slope_limit.F,v 1.6 2001/08/15 14:15:07 heimbach Exp $
C $Name:  $

#include "GMREDI_OPTIONS.h"

CStartOfInterface
      SUBROUTINE GMREDI_SLOPE_LIMIT(
     I             dSigmadRReal,
     I             depthZ,
     U             SlopeX, SlopeY,
     O             dRdSigmaLtd,
     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     |            SlopeX/Y     contains X/Y slopes              |
C     |            dRdSigmaLtd  conatins the effective dSig/dz   |
C     \==========================================================/
      IMPLICIT NONE

C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "GMREDI.h"
#include "PARAMS.h"

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 dRdSigmaLtd(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,stratLimit,f1,Smod,f2,Rnondim,Cspd,Lrho
      _RL dSigmaDrLtd
      _RL fpi
      PARAMETER(fpi=3.141592653589793047592d0)
      INTEGER i,j
      _RL     sx, sy

      if (GM_taper_scheme.EQ.'orig') then

C-      Original implementation in mitgcmuv
C       (this turns out to be the same as Cox slope clipping)
        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1

          gradSmod=SlopeX(i,j)*SlopeX(i,j)
     &            +SlopeY(i,j)*SlopeY(i,j)
          if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)

          stratLimit=-Small_Number-gradSmod*GM_rMaxSlope
          if (dSigmaDrReal(i,j).LT.stratLimit) then
            dSigmaDrLtd=dSigmaDrReal(i,j)
          else
            dSigmaDrLtd=stratLimit
          endif
          dRdSigmaLtd(i,j)=1./dSigmaDrLtd
          SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)
          SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)

         ENDDO
        ENDDO

      elseif (GM_taper_scheme.EQ.'clipping') then

C-      Cox 1987 "Slope clipping"
        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1

          gradSmod=SlopeX(i,j)*SlopeX(i,j)
     &            +SlopeY(i,j)*SlopeY(i,j)
          if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)

          dSigmaDrLtd=-(Small_Number+gradSmod*GM_rMaxSlope)
          if (dSigmaDrReal(i,j).LT.dSigmaDrLtd)
     &            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)

         ENDDO
        ENDDO

      elseif (GM_taper_scheme.EQ.'gkw91') then

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

          gradSmod=SlopeX(i,j)*SlopeX(i,j)
     &            +SlopeY(i,j)*SlopeY(i,j)
          if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)

          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

         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

          dSigmaDrLtd=dSigmaDrReal(i,j)
          if (dSigmaDrLtd.NE.0.) then
           dRdSigmaLtd(i,j)=1./dSigmaDrLtd
          else
           dRdSigmaLtd(i,j)=0.
          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

         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

          dSigmaDrLtd=dSigmaDrReal(i,j)
          if (dSigmaDrLtd.NE.0.) then
           dRdSigmaLtd(i,j)=1./dSigmaDrLtd
          else
           dRdSigmaLtd(i,j)=0.
          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*f2
          SlopeY(i,j) = sy*f1*f2

         ENDDO
        ENDDO

      elseif (GM_taper_scheme.EQ.' ') then

C-      No tapering/clipping selected so ...
        DO j=1-Oly+1,sNy+Oly-1
         DO i=1-Olx+1,sNx+Olx-1

          dSigmaDrLtd=dSigmaDrReal(i,j)
          if (dSigmaDrLtd.NE.0.) then
            dRdSigmaLtd(i,j)=1./dSigmaDrLtd
          else
            dRdSigmaLtd(i,j)=0.
          endif
          SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)
          SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)

         ENDDO
        ENDDO

      endif


#endif /* ALLOW_GMREDI */

      RETURN
      END
1	heimbach	1.7	C $Header: /u/gcmpack/models/MITgcmUV/pkg/gmredi/gmredi_slope_limit.F,v 1.6 2001/08/15 14:15:07 heimbach Exp $
2			C $Name: $
3	adcroft	1.1
4			#include "GMREDI_OPTIONS.h"
5
6			CStartOfInterface
7			SUBROUTINE GMREDI_SLOPE_LIMIT(
8			I dSigmadRReal,
9			I depthZ,
10			U SlopeX, SlopeY,
11			O dRdSigmaLtd,
12			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			C \| dSigmadRReal conatins the d/dz Sigma \|
19			C \| SlopeX/Y contains X/Y gradients of sigma \|
20			C \| depthZ conatins the height (m) of level \|
21			C \| On exit: \|
22			C \| SlopeX/Y contains X/Y slopes \|
23			C \| dRdSigmaLtd conatins the effective dSig/dz \|
24			C \==========================================================/
25			IMPLICIT NONE
26
27			C == Global variables ==
28			#include "SIZE.h"
29			#include "EEPARAMS.h"
30			#include "GMREDI.h"
31			#include "PARAMS.h"
32
33			C == Routine arguments ==
34			C
35			_RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
36			_RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37			_RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
38			_RL dRdSigmaLtd(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
39			_RL depthZ
40			INTEGER bi,bj,myThid
41			CEndOfInterface
42
43			#ifdef ALLOW_GMREDI
44
45			C == Local variables ==
46			_RL Small_Number
47	heimbach	1.7	PARAMETER(Small_Number=1.D-12)
48	adcroft	1.1	_RL gradSmod,stratLimit,f1,Smod,f2,Rnondim,Cspd,Lrho
49			_RL dSigmaDrLtd
50			_RL fpi
51	heimbach	1.6	PARAMETER(fpi=3.141592653589793047592d0)
52	adcroft	1.1	INTEGER i,j
53	heimbach	1.2	_RL sx, sy
54	adcroft	1.1
55			if (GM_taper_scheme.EQ.'orig') then
56
57			C- Original implementation in mitgcmuv
58			C (this turns out to be the same as Cox slope clipping)
59			DO j=1-Oly+1,sNy+Oly-1
60			DO i=1-Olx+1,sNx+Olx-1
61
62			gradSmod=SlopeX(i,j)*SlopeX(i,j)
63			& +SlopeY(i,j)*SlopeY(i,j)
64			if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)
65
66			stratLimit=-Small_Number-gradSmod*GM_rMaxSlope
67			if (dSigmaDrReal(i,j).LT.stratLimit) then
68			dSigmaDrLtd=dSigmaDrReal(i,j)
69			else
70			dSigmaDrLtd=stratLimit
71			endif
72			dRdSigmaLtd(i,j)=1./dSigmaDrLtd
73			SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)
74			SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)
75
76			ENDDO
77			ENDDO
78
79			elseif (GM_taper_scheme.EQ.'clipping') then
80
81			C- Cox 1987 "Slope clipping"
82			DO j=1-Oly+1,sNy+Oly-1
83			DO i=1-Olx+1,sNx+Olx-1
84
85			gradSmod=SlopeX(i,j)*SlopeX(i,j)
86			& +SlopeY(i,j)*SlopeY(i,j)
87			if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)
88
89			dSigmaDrLtd=-(Small_Number+gradSmod*GM_rMaxSlope)
90			if (dSigmaDrReal(i,j).LT.dSigmaDrLtd)
91			& dSigmaDrLtd=dSigmaDrReal(i,j)
92			dRdSigmaLtd(i,j)=1./dSigmaDrLtd
93			SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)
94			SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)
95
96			ENDDO
97			ENDDO
98
99			elseif (GM_taper_scheme.EQ.'gkw91') then
100
101			C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
102			DO j=1-Oly+1,sNy+Oly-1
103			DO i=1-Olx+1,sNx+Olx-1
104
105			gradSmod=SlopeX(i,j)*SlopeX(i,j)
106			& +SlopeY(i,j)*SlopeY(i,j)
107			if (gradSmod .NE. 0.) gradSmod=sqrt(gradSmod)
108
109			dSigmaDrLtd=dSigmaDrReal(i,j)
110			if (dSigmaDrLtd.NE.0.) then
111			dRdSigmaLtd(i,j)=1./dSigmaDrLtd
112			else
113			dRdSigmaLtd(i,j)=0.
114			endif
115			if (gradSmod.LE.GM_maxSlope*abs(dSigmaDrReal(i,j))) then
116			C If the slope is < Smax then calculate slopes properly
117			SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)
118			SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)
119			else
120			C If the slope is > Smax then taper slopes
121			SlopeX(i,j)=-SlopeX(i,j)*dSigmaDrLtd
122			& ((GM_maxSlope/gradSmod)*2)
123			SlopeY(i,j)=-SlopeY(i,j)*dSigmaDrLtd
124			& ((GM_maxSlope/gradSmod)*2)
125			endif
126
127			ENDDO
128			ENDDO
129
130			elseif (GM_taper_scheme.EQ.'dm95') then
131
132			C- Danabasoglu and McWilliams, J. Clim. 1995
133			DO j=1-Oly+1,sNy+Oly-1
134			DO i=1-Olx+1,sNx+Olx-1
135
136			dSigmaDrLtd=dSigmaDrReal(i,j)
137			if (dSigmaDrLtd.NE.0.) then
138			dRdSigmaLtd(i,j)=1./dSigmaDrLtd
139			else
140			dRdSigmaLtd(i,j)=0.
141			endif
142	heimbach	1.2	sx = -SlopeX(i,j)*dRdSigmaLtd(i,j)
143			sy = -SlopeY(i,j)*dRdSigmaLtd(i,j)
144			Smod = sxsx + sysy
145	adcroft	1.1	if (Smod.NE.0.) Smod=sqrt(Smod)
146			f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
147	heimbach	1.2	SlopeX(i,j) = sx*f1
148			SlopeY(i,j) = sy*f1
149	adcroft	1.1
150			ENDDO
151			ENDDO
152
153			elseif (GM_taper_scheme.EQ.'ldd97') then
154
155			C- Large, Danabasoglu and Doney, JPO 1997
156			DO j=1-Oly+1,sNy+Oly-1
157			DO i=1-Olx+1,sNx+Olx-1
158
159			dSigmaDrLtd=dSigmaDrReal(i,j)
160			if (dSigmaDrLtd.NE.0.) then
161			dRdSigmaLtd(i,j)=1./dSigmaDrLtd
162			else
163			dRdSigmaLtd(i,j)=0.
164			endif
165	heimbach	1.2	sx = -SlopeX(i,j)*dRdSigmaLtd(i,j)
166	adcroft	1.4	sy = -SlopeY(i,j)*dRdSigmaLtd(i,j)
167	heimbach	1.2	Smod = sxsx + sysy
168	adcroft	1.1	if (Smod.NE.0.) Smod=sqrt(Smod)
169			f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
170			Cspd=2.
171			Lrho=100.e3
172			if (FCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(Fcori(i,j,bi,bj))
173			Lrho=min(Lrho,100.e3)
174			Lrho=max(Lrho,15.e3)
175			if (Smod.NE.0.) then
176			Rnondim=depthZ/(Lrho*Smod)
177			else
178			Rnondim=0.
179			endif
180			f2=0.5(1.+sin( fpi(Rnondim-0.5)))
181	heimbach	1.5	SlopeX(i,j) = sxf1f2
182			SlopeY(i,j) = syf1f2
183	adcroft	1.1
184			ENDDO
185			ENDDO
186
187			elseif (GM_taper_scheme.EQ.' ') then
188
189			C- No tapering/clipping selected so ...
190			DO j=1-Oly+1,sNy+Oly-1
191			DO i=1-Olx+1,sNx+Olx-1
192
193			dSigmaDrLtd=dSigmaDrReal(i,j)
194			if (dSigmaDrLtd.NE.0.) then
195			dRdSigmaLtd(i,j)=1./dSigmaDrLtd
196			else
197			dRdSigmaLtd(i,j)=0.
198			endif
199			SlopeX(i,j)=-SlopeX(i,j)*dRdSigmaLtd(i,j)
200			SlopeY(i,j)=-SlopeY(i,j)*dRdSigmaLtd(i,j)
201
202			ENDDO
203			ENDDO
204
205			endif
206
207
208			#endif /* ALLOW_GMREDI */
209
210			RETURN
211			END