/[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.11 - (hide 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 mlosch 1.11 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_calc_tensor.F,v 1.10 2002/03/24 02:33:16 heimbach Exp $
2 jmc 1.8 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 heimbach 1.10 #ifdef ALLOW_AUTODIFF_TAMC
28     #include "tamc.h"
29     #include "tamc_keys.h"
30     #endif /* ALLOW_AUTODIFF_TAMC */
31    
32 adcroft 1.1 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 jmc 1.9 INTEGER bi,bj,iMin,iMax,jMin,jMax
38 adcroft 1.1 INTEGER myThid
39     CEndOfInterface
40    
41     #ifdef ALLOW_GMREDI
42    
43     C == Local variables ==
44 jmc 1.9 INTEGER i,j,k,km1,kp1
45 adcroft 1.1 _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 jmc 1.8 _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49     _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50 heimbach 1.10 _RL maskp1, maskm1, Kgm_tmp
51 adcroft 1.1
52     #ifdef GM_VISBECK_VARIABLE_K
53     _RS deltaH,zero_rs
54     PARAMETER(zero_rs=0.)
55     _RL N2,SN
56 heimbach 1.10 _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57 adcroft 1.1 #endif
58    
59 heimbach 1.10 #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 jmc 1.9 DO k=2,Nr
73     C-- 1rst loop on k : compute Tensor Coeff. at W points.
74 heimbach 1.10 km1 = MAX(1,k-1)
75     maskm1 = 1. _d 0
76     IF (k.LE.1) maskm1 = 0. _d 0
77 adcroft 1.1
78     #ifdef ALLOW_AUTODIFF_TAMC
79 heimbach 1.10 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 adcroft 1.1 #endif
93 heimbach 1.10
94 adcroft 1.1 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 heimbach 1.10 SlopeX(i,j)=0.25*( sigmaX(i+1, j ,km1) +sigmaX(i,j,km1)
99 adcroft 1.1 & +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )
100 heimbach 1.10 & *maskC(i,j,k,bi,bj)*maskm1
101     SlopeY(i,j)=0.25*( sigmaY( i ,j+1,km1) +sigmaY(i,j,km1)
102 adcroft 1.1 & +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )
103 heimbach 1.10 & *maskC(i,j,k,bi,bj)*maskm1
104     dSigmaDrReal(i,j)=sigmaR(i,j,k)*maskm1
105 adcroft 1.1
106     ENDDO
107     ENDDO
108    
109 heimbach 1.10 #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 adcroft 1.1 C Calculate slopes for use in tensor, taper and/or clip
116     CALL GMREDI_SLOPE_LIMIT(
117 jmc 1.9 U dSigmadRReal,
118 adcroft 1.1 I rF(K),
119     U SlopeX, SlopeY,
120 jmc 1.8 O SlopeSqr, taperFct,
121 adcroft 1.1 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 heimbach 1.10 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 adcroft 1.1
143 jmc 1.9 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 jmc 1.8 Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
147 adcroft 1.1
148     #ifdef GM_VISBECK_VARIABLE_K
149 jmc 1.8
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 heimbach 1.10 Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
153 jmc 1.8
154 adcroft 1.1 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 jmc 1.8 IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
168 heimbach 1.10 IF (Ssq(i,j).NE.0.) THEN
169 mlosch 1.11 N2= -Gravity*recip_RhoConst*dSigmaDrReal(i,j)
170 heimbach 1.10 SN=sqrt(Ssq(i,j)*N2)
171 heimbach 1.3 VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
172 adcroft 1.1 & *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
173 jmc 1.8 ENDIF
174 adcroft 1.1
175 jmc 1.9 #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 adcroft 1.1
184 jmc 1.9 #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 adcroft 1.1 #endif /* GM_VISBECK_VARIABLE_K */
199    
200    
201 jmc 1.9 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
202 heimbach 1.10
203 jmc 1.9 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 adcroft 1.1 #ifdef GM_VISBECK_VARIABLE_K
220 jmc 1.9 & + VisbeckK(i,j,bi,bj)
221 adcroft 1.1 #endif
222 jmc 1.9 & )*Kwz(i,j,k,bi,bj)
223 adcroft 1.1 ENDDO
224     ENDDO
225 adcroft 1.4
226 jmc 1.9 #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
227 adcroft 1.1
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 jmc 1.9 SlopeX(i,j)=sigmaX(i,j,k)
232 adcroft 1.1 & *_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 jmc 1.9 & +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )
238 adcroft 1.1 & *_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 jmc 1.9 U dSigmadRReal,
245 adcroft 1.1 I rF(K),
246     U SlopeX, SlopeY,
247 jmc 1.8 O SlopeSqr, taperFct,
248 adcroft 1.1 I bi, bj, myThid )
249    
250 jmc 1.9 #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 heimbach 1.10 & )
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 jmc 1.9 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 adcroft 1.1
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 jmc 1.9 SlopeY(i,j)=sigmaY(i,j,k)
291 adcroft 1.1 & *_maskS(i,j,k,bi,bj)
292     dSigmaDrReal(i,j)=0.25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )
293 jmc 1.9 & +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )
294 adcroft 1.1 & *_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 jmc 1.9 U dSigmadRReal,
301 adcroft 1.1 I rF(K),
302     U SlopeX, SlopeY,
303 jmc 1.8 O SlopeSqr, taperFct,
304 adcroft 1.1 I bi, bj, myThid )
305    
306 jmc 1.9 #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 heimbach 1.10 & )
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 jmc 1.9 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 adcroft 1.1 DO j=1-Oly+1,sNy+Oly-1
345     DO i=1-Olx+1,sNx+Olx-1
346 jmc 1.9 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 adcroft 1.1 ENDDO
353     ENDDO
354 jmc 1.9 GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
355     #endif /* ALLOW_TIMEAVE */
356 adcroft 1.1
357 jmc 1.9 C-- end 2nd loop on vertical level index k
358     ENDDO
359 adcroft 1.1
360    
361 jmc 1.9 #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 adcroft 1.1
370     #endif /* ALLOW_GMREDI */
371    
372     RETURN
373     END
374 heimbach 1.2
375    
376     SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
377 jmc 1.9 I bi, bj, iMin, iMax, jMin, jMax,
378 heimbach 1.2 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 jmc 1.9 INTEGER bi,bj,iMin,iMax,jMin,jMax
401 heimbach 1.2 INTEGER myThid
402     CEndOfInterface
403    
404 jmc 1.9 INTEGER i, j, k
405 heimbach 1.2
406     #ifdef ALLOW_GMREDI
407    
408 jmc 1.9 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 heimbach 1.2 ENDDO
416     ENDDO
417     #endif /* ALLOW_GMREDI */
418    
419 jmc 1.9 RETURN
420     END

  ViewVC Help
Powered by ViewVC 1.1.22