/[MITgcm]/MITgcm/pkg/gmredi/gmredi_slope_limit.F
ViewVC logotype

Contents of /MITgcm/pkg/gmredi/gmredi_slope_limit.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.6 - (show annotations) (download)
Wed Aug 15 14:15:07 2001 UTC (22 years, 10 months ago) by heimbach
Branch: MAIN
Changes since 1.5: +3 -3 lines
Added d0 to parameter fpi to make it consistent with _RL.

1 C $Header: /u/gcmpack/models/MITgcmUV/pkg/gmredi/gmredi_slope_limit.F,v 1.5 2001/08/13 18:09:23 heimbach Exp $
2 C $Name: checkpoint40pre7 $
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.141592653589793047592d0)
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 = sx*sx + sy*sy
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 sy = -SlopeY(i,j)*dRdSigmaLtd(i,j)
167 Smod = sx*sx + sy*sy
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*f2
182 SlopeY(i,j) = sy*f1*f2
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

  ViewVC Help
Powered by ViewVC 1.1.22