pkg/gmredi/gmredi_slope_limit.F

C $Header: /u/gcmpack/models/MITgcmUV/pkg/gmredi/gmredi_slope_limit.F,v 1.2 2001/01/08 20:11:04 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.E-12)
      _RL gradSmod,stratLimit,f1,Smod,f2,Rnondim,Cspd,Lrho
      _RL dSigmaDrLtd
      _RL fpi
      PARAMETER(fpi=3.141592653589793047592)
      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)
          sx   = -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
          SlopeY(i,j) = sy*f1

         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	C $Header: /u/gcmpack/models/MITgcmUV/pkg/gmredi/gmredi_slope_limit.F,v 1.2 2001/01/08 20:11:04 heimbach Exp $
2	C $Name: $
3
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	PARAMETER(Small_Number=1.E-12)
48	_RL gradSmod,stratLimit,f1,Smod,f2,Rnondim,Cspd,Lrho
49	_RL dSigmaDrLtd
50	_RL fpi
51	PARAMETER(fpi=3.141592653589793047592)
52	INTEGER i,j
53	_RL sx, sy
54
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	sx = -SlopeX(i,j)*dRdSigmaLtd(i,j)
143	sy = -SlopeY(i,j)*dRdSigmaLtd(i,j)
144	Smod = sxsx + sysy
145	if (Smod.NE.0.) Smod=sqrt(Smod)
146	f1=0.5*(1.+tanh( (GM_Scrit-Smod)/GM_Sd ))
147	SlopeX(i,j) = sx*f1
148	SlopeY(i,j) = sy*f1
149
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	sx = -SlopeX(i,j)*dRdSigmaLtd(i,j)
166	sx = -SlopeY(i,j)*dRdSigmaLtd(i,j)
167	Smod = sxsx + sysy
168	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	SlopeX(i,j) = sx*f1
182	SlopeY(i,j) = sy*f1
183
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