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

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

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


Revision 1.11 - (show annotations) (download)
Wed Sep 25 19:36:50 2002 UTC (21 years, 8 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint46n_post, checkpoint46l_post, checkpoint46l_pre, checkpoint46j_pre, checkpoint46j_post, checkpoint46k_post, checkpoint46m_post, checkpoint46i_post, checkpoint46h_post
Changes since 1.10: +2 -2 lines
o cleaned up the use of rhoNil and rhoConst.
  - rhoNil should only appear in the LINEAR equation of state, everywhere
    else rhoNil is replaced by rhoConst, e.g. find_rho computes rho-rhoConst
    and the dynamical equations are all divided by rhoConst
o introduced new parameter rhoConstFresh, a reference density of fresh
  water, to remove the fresh water flux's dependence on rhoNil. The default
  value is 999.8 kg/m^3
o cleanup up external_forcing.F and external_forcing_surf.F
  - can now be used by both OCEANIC and OCEANICP

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

  ViewVC Help
Powered by ViewVC 1.1.22