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

Annotation of /MITgcm/pkg/gmredi/gmredi_slope_psi.F

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

Revision 1.6 - (hide annotations) (download)
Tue Jan 21 19:34:13 2003 UTC (21 years, 3 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint51k_post, checkpoint52l_pre, hrcube4, hrcube5, checkpoint50c_post, checkpoint52d_pre, checkpoint48e_post, checkpoint50c_pre, checkpoint52j_pre, checkpoint51o_pre, checkpoint54d_post, checkpoint54e_post, checkpoint51l_post, checkpoint48i_post, checkpoint52l_post, checkpoint52k_post, checkpoint55, checkpoint54, checkpoint56, checkpoint51, checkpoint50, checkpoint53, checkpoint52, checkpoint50d_post, checkpoint52f_post, checkpoint50b_pre, checkpoint54f_post, checkpoint51f_post, checkpoint48b_post, checkpoint51d_post, checkpoint48c_pre, checkpoint51t_post, checkpoint51n_post, checkpoint55i_post, checkpoint52i_pre, hrcube_1, hrcube_2, hrcube_3, checkpoint51s_post, checkpoint55c_post, checkpoint48d_pre, checkpoint51j_post, checkpoint52e_pre, checkpoint52e_post, checkpoint51n_pre, checkpoint53d_post, checkpoint48d_post, checkpoint48f_post, checkpoint52b_pre, checkpoint54b_post, checkpoint51l_pre, checkpoint52m_post, checkpoint55g_post, checkpoint48h_post, checkpoint51q_post, checkpoint51b_pre, checkpoint52b_post, checkpoint52c_post, checkpoint51h_pre, checkpoint48a_post, checkpoint50f_post, checkpoint50a_post, checkpoint50f_pre, checkpoint52f_pre, checkpoint55d_post, checkpoint47j_post, checkpoint54a_pre, checkpoint53c_post, checkpoint55d_pre, checkpoint55j_post, branchpoint-genmake2, checkpoint54a_post, checkpoint55h_post, checkpoint51r_post, checkpoint48c_post, checkpoint51i_post, checkpoint55b_post, checkpoint51b_post, checkpoint51c_post, checkpoint53a_post, checkpoint55f_post, checkpoint52d_post, checkpoint53g_post, checkpoint50g_post, checkpoint52a_pre, checkpoint50h_post, checkpoint52i_post, checkpoint50e_pre, checkpoint50i_post, checkpoint51i_pre, checkpoint52h_pre, checkpoint53f_post, checkpoint52j_post, checkpoint50e_post, branch-netcdf, checkpoint50d_pre, checkpoint52n_post, checkpoint53b_pre, checkpoint51e_post, checkpoint55a_post, checkpoint48, checkpoint49, checkpoint51o_post, checkpoint51f_pre, checkpoint48g_post, checkpoint53b_post, checkpoint52a_post, checkpoint51g_post, ecco_c52_e35, checkpoint50b_post, checkpoint51m_post, checkpoint53d_pre, checkpoint55e_post, checkpoint54c_post, checkpoint51a_post, checkpoint51p_post, checkpoint51u_post
Branch point for: branch-genmake2, branch-nonh, tg2-branch, netcdf-sm0, checkpoint51n_branch
Changes since 1.5: +69 -19 lines 
Yet more changes:
o adgmredi_calc_tensor
  avoiding all recomputation of gmredi_slope_limit
o adgmredi_x/y/rtransport
  added flag for excessive storing to avoid recomp. of
  u/v/rtans, dTdx/y/z
  -> this is not really necessary and very memory-consuming
o adgmredi_slope_psi:
  consistency with gmredi_slope_limit in treatment of GM_slopeSqCutoff
o gmredi_slope_limit
  re-activated full calculation of taperfct for case 'ac02'
1 heimbach 1.6 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_slope_psi.F,v 1.5 2003/01/12 21:27:20 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "GMREDI_OPTIONS.h"
5    
6     CStartOfInterface
7 jmc 1.5 SUBROUTINE GMREDI_SLOPE_PSI(
8 jmc 1.1 I dSigmaDrW,dSigmaDrS,
9 heimbach 1.3 I depthZ,K,
10 jmc 1.1 U SlopeX, SlopeY,
11     O taperX, taperY,
12     I bi,bj, myThid )
13     C /==========================================================\
14 jmc 1.5 C | SUBROUTINE GMREDI_SLOPE_PSI |
15 jmc 1.1 C | o Calculate slopes for use in GM/Redi tensor |
16     C |==========================================================|
17     C | On entry: |
18 heimbach 1.3 C | dSigmaDrW conatins the d/dz Sigma |
19 jmc 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 heimbach 1.3 C | dSigmaDrW conatins the effective dSig/dz |
23 jmc 1.1 C | SlopeX/Y contains X/Y slopes |
24     C | taperFct contains tapering funct. value ; |
25     C | = 1 when using no tapering |
26     C \==========================================================/
27     IMPLICIT NONE
28    
29     C == Global variables ==
30     #include "SIZE.h"
31     #include "EEPARAMS.h"
32     #include "GMREDI.h"
33     #include "PARAMS.h"
34    
35 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
36     #include "tamc.h"
37     #include "tamc_keys.h"
38     #endif /* ALLOW_AUTODIFF_TAMC */
39    
40 jmc 1.1 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 dSigmaDrW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45     _RL dSigmaDrS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46     _RL taperX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47     _RL taperY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48     _RL depthZ
49 heimbach 1.3 INTEGER K,bi,bj,myThid
50 jmc 1.1 CEndOfInterface
51    
52     #ifdef ALLOW_GMREDI
53 heimbach 1.3 #ifdef GM_BOLUS_ADVEC
54 jmc 1.1
55     C == Local variables ==
56 heimbach 1.3 _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 heimbach 1.6 _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60 heimbach 1.3 _RL f1,Smod,f2,Rnondim,Cspd,Lrho
61 jmc 1.1 _RL maxSlopeSqr
62 heimbach 1.6 _RL tmpvar
63 jmc 1.1 _RL fpi
64     PARAMETER(fpi=3.141592653589793047592d0)
65     INTEGER i,j
66    
67 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
68     act1 = bi - myBxLo(myThid)
69     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
70     act2 = bj - myByLo(myThid)
71     max2 = myByHi(myThid) - myByLo(myThid) + 1
72     act3 = myThid - 1
73     max3 = nTx*nTy
74     act4 = ikey_dynamics - 1
75 heimbach 1.6 igmkey = (act1 + 1) + act2*max1
76     & + act3*max1*max2
77     & + act4*max1*max2*max3
78     kkey = (igmkey-1)*Nr + k
79 heimbach 1.3 #endif /* ALLOW_AUTODIFF_TAMC */
80    
81 jmc 1.1 IF (GM_taper_scheme.EQ.'orig' .OR.
82     & GM_taper_scheme.EQ.'clipping') THEN
83    
84 heimbach 1.6 #ifdef GM_EXCLUDE_CLIPPING
85    
86     STOP 'Need to compile without "#define GM_EXCLUDE_CLIPPING"'
87    
88     #else /* GM_EXCLUDE_CLIPPING */
89    
90 jmc 1.1 C- Original implementation in mitgcmuv
91     C (this turns out to be the same as Cox slope clipping)
92    
93 heimbach 1.3 C-- X-comp
94    
95     #ifdef ALLOW_AUTODIFF_TAMC
96     DO j=1-Oly+1,sNy+Oly-1
97     DO i=1-Olx+1,sNx+Olx-1
98     dSigmaDrLtd(i,j) = 0. _d 0
99     ENDDO
100     ENDDO
101     #endif /* ALLOW_AUTODIFF_TAMC */
102    
103 jmc 1.1 C- Cox 1987 "Slope clipping"
104     DO j=1-Oly+1,sNy+Oly-1
105     DO i=1-Olx+1,sNx+Olx-1
106 heimbach 1.3 dsigmadrltd(i,j) = -(GM_Small_Number+
107     & abs(SlopeX(i,j))*GM_rMaxSlope)
108     ENDDO
109     ENDDO
110     #ifdef ALLOW_AUTODIFF_TAMC
111     CADJ STORE dSigmaDrltd(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
112     CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
113     #endif
114     DO j=1-Oly+1,sNy+Oly-1
115     DO i=1-Olx+1,sNx+Olx-1
116     IF (dSigmaDrW(i,j).GE.dsigmadrltd(i,j))
117     & dSigmaDrW(i,j) = dsigmadrltd(i,j)
118     ENDDO
119     ENDDO
120     #ifdef ALLOW_AUTODIFF_TAMC
121     CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
122     #endif
123     DO j=1-Oly+1,sNy+Oly-1
124     DO i=1-Olx+1,sNx+Olx-1
125     SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
126     taperX(i,j)=1. _d 0
127     ENDDO
128     ENDDO
129 jmc 1.1
130 heimbach 1.3 C-- Y-comp
131 jmc 1.1
132 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
133     DO j=1-Oly+1,sNy+Oly-1
134     DO i=1-Olx+1,sNx+Olx-1
135     dSigmaDrLtd(i,j) = 0. _d 0
136     ENDDO
137     ENDDO
138     #endif /* ALLOW_AUTODIFF_TAMC */
139     DO j=1-Oly+1,sNy+Oly-1
140     DO i=1-Olx+1,sNx+Olx-1
141     dsigmadrltd(i,j) = -(GM_Small_Number+
142     & abs(SlopeY(i,j))*GM_rMaxSlope)
143     ENDDO
144     ENDDO
145     #ifdef ALLOW_AUTODIFF_TAMC
146     CADJ STORE dSigmaDrltd(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
147     CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
148     #endif
149     DO j=1-Oly+1,sNy+Oly-1
150     DO i=1-Olx+1,sNx+Olx-1
151     IF (dSigmaDrS(i,j).GE.dsigmadrltd(i,j))
152     & dSigmaDrS(i,j) = dsigmadrltd(i,j)
153     ENDDO
154     ENDDO
155     #ifdef ALLOW_AUTODIFF_TAMC
156     CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
157     #endif
158     DO j=1-Oly+1,sNy+Oly-1
159     DO i=1-Olx+1,sNx+Olx-1
160 jmc 1.1 SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
161     taperY(i,j)=1. _d 0
162     ENDDO
163     ENDDO
164    
165 heimbach 1.6 #endif /* GM_EXCLUDE_CLIPPING */
166    
167 jmc 1.1 ELSE
168    
169 heimbach 1.6 #ifdef GM_EXCLUDE_TAPERING
170    
171     STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
172    
173     #else /* GM_EXCLUDE_TAPERING */
174    
175 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
176     CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
177     CADJ STORE dSigmaDrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
178     #endif
179    
180 heimbach 1.6 tmpvar = sqrt( GM_slopeSqCutoff )
181    
182 jmc 1.1 C- Compute the slope, no clipping, but avoid reverse slope in negatively
183     C stratified (Sigma_Z > 0) region :
184     DO j=1-Oly+1,sNy+Oly-1
185     DO i=1-Olx+1,sNx+Olx-1
186 heimbach 1.3 dsigmadrltd(i,j) = -GM_Small_Number
187     IF (dSigmaDrW(i,j).GE.dsigmadrltd(i,j))
188     & dSigmaDrW(i,j) = dsigmadrltd(i,j)
189     ENDDO
190     ENDDO
191     #ifdef ALLOW_AUTODIFF_TAMC
192     CADJ STORE dsigmadrW(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
193     #endif
194     DO j=1-Oly+1,sNy+Oly-1
195     DO i=1-Olx+1,sNx+Olx-1
196     drdsigmaltd(i,j) = 1./dSigmaDrW(i,j)
197     SlopeX(i,j) = -SlopeX(i,j)/dSigmaDrW(i,j)
198 heimbach 1.6 taperX(i,j)= 1. _d 0
199     ENDDO
200     ENDDO
201     #ifdef ALLOW_AUTODIFF_TAMC
202     CADJ STORE slopex(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
203     #endif
204     DO j=1-Oly+1,sNy+Oly-1
205     DO i=1-Olx+1,sNx+Olx-1
206     IF ( ABS(SlopeX(i,j)) .GE. tmpvar ) THEN
207     SlopeX(i,j) = SIGN(tmpvar,SlopeX(i,j))
208     taperX(i,j) = 0. _d 0
209     ENDIF
210 heimbach 1.3 ENDDO
211     ENDDO
212 jmc 1.1
213 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
214     CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
215     CADJ STORE dSigmaDrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
216     #endif
217 jmc 1.1
218 heimbach 1.3 DO j=1-Oly+1,sNy+Oly-1
219     DO i=1-Olx+1,sNx+Olx-1
220     dsigmadrltd(i,j) = -GM_Small_Number
221     IF (dSigmaDrS(i,j).GE.dsigmadrltd(i,j))
222     & dSigmaDrS(i,j) = dsigmadrltd(i,j)
223     ENDDO
224     ENDDO
225     #ifdef ALLOW_AUTODIFF_TAMC
226     CADJ STORE dsigmadrS(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
227     #endif
228     DO j=1-Oly+1,sNy+Oly-1
229     DO i=1-Olx+1,sNx+Olx-1
230 jmc 1.1 SlopeY(i,j) = -SlopeY(i,j)/dSigmaDrS(i,j)
231     taperY(i,j)=1. _d 0
232     ENDDO
233     ENDDO
234 heimbach 1.6 #ifdef ALLOW_AUTODIFF_TAMC
235     CADJ STORE slopey(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
236     #endif
237     DO j=1-Oly+1,sNy+Oly-1
238     DO i=1-Olx+1,sNx+Olx-1
239     IF ( ABS(SlopeY(i,j)) .GE. tmpvar ) THEN
240     SlopeY(i,j) = SIGN(tmpvar,SlopeY(i,j))
241     taperY(i,j) = 0. _d 0
242     ENDIF
243     ENDDO
244     ENDDO
245 jmc 1.1
246     C- Compute the tapering function for the GM+Redi tensor :
247    
248 heimbach 1.3 #ifdef ALLOW_AUTODIFF_TAMC
249     CADJ STORE slopeX(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
250     CADJ STORE slopeY(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
251     #endif
252    
253 jmc 1.1 IF (GM_taper_scheme.EQ.'linear') THEN
254    
255     C- Simplest adiabatic tapering = Smax/Slope (linear)
256     DO j=1-Oly+1,sNy+Oly-1
257     DO i=1-Olx+1,sNx+Olx-1
258    
259 heimbach 1.6 Smod = abs(SlopeX(i,j))
260     IF ( Smod .GT. GM_maxSlope .AND.
261     & Smod .LT. tmpvar )
262     & taperX(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
263     Smod = abs(SlopeY(i,j))
264     IF ( Smod .GT. GM_maxSlope .AND.
265     & Smod .LT. tmpvar )
266     & taperY(i,j)=GM_maxSlope/(Smod+GM_Small_Number)
267 jmc 1.1
268     ENDDO
269     ENDDO
270    
271     ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
272    
273     C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
274 heimbach 1.6
275 jmc 1.1 maxSlopeSqr = GM_maxSlope*GM_maxSlope
276     DO j=1-Oly+1,sNy+Oly-1
277     DO i=1-Olx+1,sNx+Olx-1
278    
279 heimbach 1.6 IF ( abs(SlopeX(i,j)) .GT. GM_maxSlope .AND.
280     & abs(SlopeX(i,j)) .LT. tmpvar )
281     & taperX(i,j)=maxSlopeSqr/
282     & ( SlopeX(i,j)*SlopeX(i,j) + GM_Small_Number )
283     IF ( abs(SlopeY(i,j)) .GT. GM_maxSlope .AND.
284     & abs(SlopeY(i,j)) .LT. tmpvar )
285     & taperY(i,j)=maxSlopeSqr/
286     & ( SlopeY(i,j)*SlopeY(i,j) + GM_Small_Number )
287 jmc 1.1
288     ENDDO
289     ENDDO
290    
291     ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
292    
293     C- Danabasoglu and McWilliams, J. Clim. 1995
294     DO j=1-Oly+1,sNy+Oly-1
295     DO i=1-Olx+1,sNx+Olx-1
296    
297     Smod = abs(SlopeX(i,j))
298 heimbach 1.4 taperX(i,j)=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
299 jmc 1.1 Smod = abs(SlopeY(i,j))
300 heimbach 1.4 taperY(i,j)=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
301 jmc 1.1
302     ENDDO
303     ENDDO
304    
305     ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
306    
307     C- Large, Danabasoglu and Doney, JPO 1997
308     DO j=1-Oly+1,sNy+Oly-1
309     DO i=1-Olx+1,sNx+Olx-1
310    
311 heimbach 1.4 Cspd=2. _d 0
312 jmc 1.5 Lrho=100. _d 3
313     if (fCori(i,j,bi,bj).NE.0.) Lrho=Cspd/abs(fCori(i,j,bi,bj))
314 heimbach 1.2 Lrho=min(Lrho , 100. _d 3)
315     Lrho=max(Lrho , 15. _d 3)
316 jmc 1.1
317     Smod = abs(SlopeX(i,j))
318 heimbach 1.6 if ( Smod .LT. tmpvar ) then
319 heimbach 1.4 f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
320 jmc 1.1 if (Smod.NE.0.) then
321     Rnondim=depthZ/(Lrho*Smod)
322     else
323     Rnondim=0.
324     endif
325 heimbach 1.4 f2=op5*( 1. _d 0 + sin( fpi*(Rnondim-op5)))
326 jmc 1.1 taperX(i,j)=f1*f2
327 heimbach 1.6 endif
328 jmc 1.1
329     Smod = abs(SlopeY(i,j))
330 heimbach 1.6 if ( Smod .LT. tmpvar ) then
331 jmc 1.5 f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
332 jmc 1.1 if (Smod.NE.0.) then
333     Rnondim=depthZ/(Lrho*Smod)
334     else
335     Rnondim=0.
336     endif
337 jmc 1.5 f2=op5*( 1. _d 0 + sin( fpi*(Rnondim-op5)))
338 jmc 1.1 taperY(i,j)=f1*f2
339 heimbach 1.6 endif
340 jmc 1.1
341     ENDDO
342     ENDDO
343    
344     ELSEIF (GM_taper_scheme.NE.' ') THEN
345     STOP 'GMREDI_SLOPE_PSI: Bad GM_taper_scheme'
346     ENDIF
347 heimbach 1.6
348     #endif /* GM_EXCLUDE_TAPERING */
349 jmc 1.1
350     ENDIF
351    
352    
353 heimbach 1.3 #endif /* BOLUS_ADVEC */
354 jmc 1.1 #endif /* ALLOW_GMREDI */
355    
356     RETURN
357     END
  ViewVC Help
Powered by ViewVC 1.1.22