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

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

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


Revision 1.9 - (hide annotations) (download)
Sun Dec 16 18:54:49 2001 UTC (22 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint44e_post, checkpoint44f_post, checkpoint43a-release1mods, chkpt44d_post, checkpoint44e_pre, release1-branch_tutorials, chkpt44a_post, checkpoint44h_pre, chkpt44c_pre, checkpoint44g_post, release1-branch-end, checkpoint44b_post, chkpt44a_pre, checkpoint44b_pre, checkpoint44, chkpt44c_post, checkpoint44f_pre, release1-branch_branchpoint
Branch point for: release1_final, release1-branch
Changes since 1.8: +152 -55 lines
Modification to the GMREDI package :
 change units of tensor-K arrays, scale now like diffusivity
 initialise all common block arrays in S/R gmredi_init
 add option to use different isopycnal(Redi) & GM diffusivity
 add option to use the advective GM form or the skew-flux form (=default)
 bug in non_unity_diagonal part fixed.

1 jmc 1.8 C $Header: /u/gcmpack/models/MITgcmUV/pkg/gmredi/gmredi_calc_tensor.F,v 1.7 2001/08/21 15:27:19 heimbach Exp $
2     C $Name: $
3 adcroft 1.1
4     #include "GMREDI_OPTIONS.h"
5    
6     CStartOfInterface
7     SUBROUTINE GMREDI_CALC_TENSOR(
8 jmc 1.9 I bi, bj, iMin, iMax, jMin, jMax,
9 adcroft 1.1 I sigmaX, sigmaY, sigmaR,
10     I myThid )
11     C /==========================================================\
12     C | SUBROUTINE GMREDI_CALC_TENSOR |
13     C | o Calculate tensor elements for GM/Redi tensor. |
14     C |==========================================================|
15     C \==========================================================/
16     IMPLICIT NONE
17    
18     C == Global variables ==
19     #include "SIZE.h"
20     #include "GRID.h"
21     #include "DYNVARS.h"
22     #include "EEPARAMS.h"
23     #include "PARAMS.h"
24     #include "GMREDI.h"
25     #include "GMREDI_DIAGS.h"
26    
27     C == Routine arguments ==
28     C
29     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
30     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
31     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
32 jmc 1.9 INTEGER bi,bj,iMin,iMax,jMin,jMax
33 adcroft 1.1 INTEGER myThid
34     CEndOfInterface
35    
36     #ifdef ALLOW_GMREDI
37    
38     C == Local variables ==
39 jmc 1.9 INTEGER i,j,k,km1,kp1
40 adcroft 1.1 _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41     _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
42     _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43 jmc 1.8 _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44     _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45 jmc 1.9 _RL maskp1, Kgm_tmp
46 adcroft 1.1
47     #ifdef GM_VISBECK_VARIABLE_K
48     _RS deltaH,zero_rs
49     PARAMETER(zero_rs=0.)
50     _RL N2,SN
51 jmc 1.8 _RL Ssq
52 adcroft 1.1 #endif
53    
54 jmc 1.9 DO k=2,Nr
55     C-- 1rst loop on k : compute Tensor Coeff. at W points.
56     c km1 = MAX(1,k-1)
57 adcroft 1.1
58     #ifdef ALLOW_AUTODIFF_TAMC
59     !HPF$ INDEPENDENT
60     #endif
61     DO j=1-Oly+1,sNy+Oly-1
62     #ifdef ALLOW_AUTODIFF_TAMC
63     !HPF$ INDEPENDENT
64     #endif
65     DO i=1-Olx+1,sNx+Olx-1
66    
67     C Gradient of Sigma at rVel points
68 jmc 1.9 SlopeX(i,j)=0.25*( sigmaX(i+1, j ,k-1) +sigmaX(i,j,k-1)
69 adcroft 1.1 & +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )
70 jmc 1.8 & *maskC(i,j,k,bi,bj)
71 jmc 1.9 SlopeY(i,j)=0.25*( sigmaY( i ,j+1,k-1) +sigmaY(i,j,k-1)
72 adcroft 1.1 & +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )
73 jmc 1.8 & *maskC(i,j,k,bi,bj)
74 adcroft 1.1 dSigmaDrReal(i,j)=sigmaR(i,j,k)
75    
76     ENDDO
77     ENDDO
78    
79     C Calculate slopes for use in tensor, taper and/or clip
80     CALL GMREDI_SLOPE_LIMIT(
81 jmc 1.9 U dSigmadRReal,
82 adcroft 1.1 I rF(K),
83     U SlopeX, SlopeY,
84 jmc 1.8 O SlopeSqr, taperFct,
85 adcroft 1.1 I bi, bj, myThid )
86    
87     DO j=1-Oly+1,sNy+Oly-1
88     DO i=1-Olx+1,sNx+Olx-1
89    
90     C Mask Iso-neutral slopes
91 jmc 1.8 SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
92     SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
93     SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
94     c Ssq=SlopeX(i,j)*SlopeX(i,j)+SlopeY(i,j)*SlopeY(i,j)
95 adcroft 1.1
96 jmc 1.9 C Components of Redi/GM tensor
97     Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
98     Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
99 jmc 1.8 Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
100 adcroft 1.1
101     #ifdef GM_VISBECK_VARIABLE_K
102 jmc 1.8
103     C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
104     C but don't know if *taperFct (or **2 ?) is necessary
105     Ssq=SlopeSqr(i,j)*taperFct(i,j)
106    
107 adcroft 1.1 C-- Depth average of M^2/N^2 * N
108    
109     C Calculate terms for mean Richardson number
110     C which is used in the "variable K" parameterisaton.
111     C Distance between interface above layer and the integration depth
112     deltaH=abs(GM_Visbeck_depth)-abs(rF(k))
113     C If positive we limit this to the layer thickness
114     deltaH=min(deltaH,drF(k))
115     C If negative then we are below the integration level
116     deltaH=max(deltaH,zero_rs)
117     C Now we convert deltaH to a non-dimensional fraction
118     deltaH=deltaH/GM_Visbeck_depth
119    
120 jmc 1.8 IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
121     IF (Ssq.NE.0.) THEN
122     N2= -Gravity*recip_Rhonil*dSigmaDrReal(i,j)
123 adcroft 1.1 SN=sqrt(Ssq*N2)
124 heimbach 1.3 VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
125 adcroft 1.1 & *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
126 jmc 1.8 ENDIF
127 adcroft 1.1
128 jmc 1.9 #endif /* GM_VISBECK_VARIABLE_K */
129    
130     ENDDO
131     ENDDO
132    
133     C-- end 1rst loop on vertical level index k
134     ENDDO
135    
136 adcroft 1.1
137 jmc 1.9 #ifdef GM_VISBECK_VARIABLE_K
138     IF ( GM_Visbeck_alpha.NE.0. ) THEN
139     C- Limit range that KapGM can take
140     DO j=1-Oly+1,sNy+Oly-1
141     DO i=1-Olx+1,sNx+Olx-1
142     VisbeckK(i,j,bi,bj)=
143     & MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
144     #ifdef ALLOW_TIMEAVE
145     Visbeck_K_T(i,j,bi,bj)=Visbeck_K_T(i,j,bi,bj)
146     & +VisbeckK(i,j,bi,bj)*deltaTclock
147     #endif
148     ENDDO
149     ENDDO
150     ENDIF
151 adcroft 1.1 #endif /* GM_VISBECK_VARIABLE_K */
152    
153    
154 jmc 1.9 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
155     C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.
156     DO k=1,Nr
157     kp1 = MIN(Nr,k+1)
158     maskp1 = 1. _d 0
159     IF (k.GE.Nr) maskp1 = 0. _d 0
160    
161     C- express the Tensor in term of Diffusivity (= m**2 / s )
162     DO j=1-Oly+1,sNy+Oly-1
163     DO i=1-Olx+1,sNx+Olx-1
164     Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
165     #ifdef GM_VISBECK_VARIABLE_K
166     & + VisbeckK(i,j,bi,bj)*(1.+GM_skewflx)
167     #endif
168     Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
169     Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
170     Kwz(i,j,k,bi,bj)= ( GM_isopycK
171 adcroft 1.1 #ifdef GM_VISBECK_VARIABLE_K
172 jmc 1.9 & + VisbeckK(i,j,bi,bj)
173 adcroft 1.1 #endif
174 jmc 1.9 & )*Kwz(i,j,k,bi,bj)
175 adcroft 1.1 ENDDO
176     ENDDO
177 adcroft 1.4
178 jmc 1.9 #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
179 adcroft 1.1
180     C Gradient of Sigma at U points
181     DO j=1-Oly+1,sNy+Oly-1
182     DO i=1-Olx+1,sNx+Olx-1
183 jmc 1.9 SlopeX(i,j)=sigmaX(i,j,k)
184 adcroft 1.1 & *_maskW(i,j,k,bi,bj)
185     SlopeY(i,j)=0.25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)
186     & +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )
187     & *_maskW(i,j,k,bi,bj)
188     dSigmaDrReal(i,j)=0.25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )
189 jmc 1.9 & +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )
190 adcroft 1.1 & *_maskW(i,j,k,bi,bj)
191     ENDDO
192     ENDDO
193    
194     C Calculate slopes for use in tensor, taper and/or clip
195     CALL GMREDI_SLOPE_LIMIT(
196 jmc 1.9 U dSigmadRReal,
197 adcroft 1.1 I rF(K),
198     U SlopeX, SlopeY,
199 jmc 1.8 O SlopeSqr, taperFct,
200 adcroft 1.1 I bi, bj, myThid )
201    
202 jmc 1.9 #ifdef GM_NON_UNITY_DIAGONAL
203     DO j=1-Oly+1,sNy+Oly-1
204     DO i=1-Olx+1,sNx+Olx-1
205     Kux(i,j,k,bi,bj) =
206     & ( GM_isopycK
207     #ifdef GM_VISBECK_VARIABLE_K
208     & +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
209     #endif
210     & )*taperFct(i,j)
211     Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
212     ENDDO
213     ENDDO
214     #endif /* GM_NON_UNITY_DIAGONAL */
215    
216     #ifdef GM_EXTRA_DIAGONAL
217     IF (GM_ExtraDiag) THEN
218     DO j=1-Oly+1,sNy+Oly-1
219     DO i=1-Olx+1,sNx+Olx-1
220     Kuz(i,j,k,bi,bj) =
221     & ( GM_isopycK - GM_skewflx*GM_background_K
222     #ifdef GM_VISBECK_VARIABLE_K
223     & +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
224     #endif
225     & )*SlopeX(i,j)*taperFct(i,j)
226     ENDDO
227     ENDDO
228     ENDIF
229     #endif /* GM_EXTRA_DIAGONAL */
230 adcroft 1.1
231     C Gradient of Sigma at V points
232     DO j=1-Oly+1,sNy+Oly-1
233     DO i=1-Olx+1,sNx+Olx-1
234     SlopeX(i,j)=0.25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
235     & +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )
236     & *_maskS(i,j,k,bi,bj)
237 jmc 1.9 SlopeY(i,j)=sigmaY(i,j,k)
238 adcroft 1.1 & *_maskS(i,j,k,bi,bj)
239     dSigmaDrReal(i,j)=0.25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )
240 jmc 1.9 & +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )
241 adcroft 1.1 & *_maskS(i,j,k,bi,bj)
242     ENDDO
243     ENDDO
244    
245     C Calculate slopes for use in tensor, taper and/or clip
246     CALL GMREDI_SLOPE_LIMIT(
247 jmc 1.9 U dSigmadRReal,
248 adcroft 1.1 I rF(K),
249     U SlopeX, SlopeY,
250 jmc 1.8 O SlopeSqr, taperFct,
251 adcroft 1.1 I bi, bj, myThid )
252    
253 jmc 1.9 #ifdef GM_NON_UNITY_DIAGONAL
254     DO j=1-Oly+1,sNy+Oly-1
255     DO i=1-Olx+1,sNx+Olx-1
256     Kvy(i,j,k,bi,bj) =
257     & ( GM_isopycK
258     #ifdef GM_VISBECK_VARIABLE_K
259     & +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
260     #endif
261     & )*taperFct(i,j)
262     Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
263     ENDDO
264     ENDDO
265     #endif /* GM_NON_UNITY_DIAGONAL */
266    
267     #ifdef GM_EXTRA_DIAGONAL
268     IF (GM_ExtraDiag) THEN
269     DO j=1-Oly+1,sNy+Oly-1
270     DO i=1-Olx+1,sNx+Olx-1
271     Kvz(i,j,k,bi,bj) =
272     & ( GM_isopycK - GM_skewflx*GM_background_K
273     #ifdef GM_VISBECK_VARIABLE_K
274     & +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
275     #endif
276     & )*SlopeY(i,j)*taperFct(i,j)
277     ENDDO
278     ENDDO
279     ENDIF
280     #endif /* GM_EXTRA_DIAGONAL */
281    
282     #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
283    
284     #ifdef ALLOW_TIMEAVE
285     C-- Time-average
286 adcroft 1.1 DO j=1-Oly+1,sNy+Oly-1
287     DO i=1-Olx+1,sNx+Olx-1
288 jmc 1.9 GM_Kwx_T(i,j,k,bi,bj)=GM_Kwx_T(i,j,k,bi,bj)
289     & +Kwx(i,j,k,bi,bj)*deltaTclock
290     GM_Kwy_T(i,j,k,bi,bj)=GM_Kwy_T(i,j,k,bi,bj)
291     & +Kwy(i,j,k,bi,bj)*deltaTclock
292     GM_Kwz_T(i,j,k,bi,bj)=GM_Kwz_T(i,j,k,bi,bj)
293     & +Kwz(i,j,k,bi,bj)*deltaTclock
294 adcroft 1.1 ENDDO
295     ENDDO
296 jmc 1.9 GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
297     #endif /* ALLOW_TIMEAVE */
298 adcroft 1.1
299 jmc 1.9 C-- end 2nd loop on vertical level index k
300     ENDDO
301 adcroft 1.1
302    
303 jmc 1.9 #ifdef GM_BOLUS_ADVEC
304     IF (GM_AdvForm) THEN
305     CALL GMREDI_CALC_PSI_B(
306     I bi, bj, iMin, iMax, jMin, jMax,
307     I sigmaX, sigmaY, sigmaR,
308     I myThid )
309     ENDIF
310     #endif
311 adcroft 1.1
312     #endif /* ALLOW_GMREDI */
313    
314     RETURN
315     END
316 heimbach 1.2
317    
318     SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
319 jmc 1.9 I bi, bj, iMin, iMax, jMin, jMax,
320 heimbach 1.2 I sigmaX, sigmaY, sigmaR,
321     I myThid )
322     C /==========================================================\
323     C | SUBROUTINE GMREDI_CALC_TENSOR |
324     C | o Calculate tensor elements for GM/Redi tensor. |
325     C |==========================================================|
326     C \==========================================================/
327     IMPLICIT NONE
328    
329     C == Global variables ==
330     #include "SIZE.h"
331     #include "GRID.h"
332     #include "DYNVARS.h"
333     #include "EEPARAMS.h"
334     #include "PARAMS.h"
335     #include "GMREDI.h"
336    
337     C == Routine arguments ==
338     C
339     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
340     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
341     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
342 jmc 1.9 INTEGER bi,bj,iMin,iMax,jMin,jMax
343 heimbach 1.2 INTEGER myThid
344     CEndOfInterface
345    
346 jmc 1.9 INTEGER i, j, k
347 heimbach 1.2
348     #ifdef ALLOW_GMREDI
349    
350 jmc 1.9 DO k=1,Nr
351     DO j=1-Oly+1,sNy+Oly-1
352     DO i=1-Olx+1,sNx+Olx-1
353     Kwx(i,j,k,bi,bj) = 0.0
354     Kwy(i,j,k,bi,bj) = 0.0
355     Kwz(i,j,k,bi,bj) = 0.0
356     ENDDO
357 heimbach 1.2 ENDDO
358     ENDDO
359     #endif /* ALLOW_GMREDI */
360    
361 jmc 1.9 RETURN
362     END

  ViewVC Help
Powered by ViewVC 1.1.22