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

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

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

revision 1.8 by jmc, Tue Dec 4 15:09:34 2001 UTC revision 1.9 by jmc, Sun Dec 16 18:54:49 2001 UTC
# Line 5  C $Name$ Line 5  C $Name$
5    
6  CStartOfInterface  CStartOfInterface
7        SUBROUTINE GMREDI_CALC_TENSOR(        SUBROUTINE GMREDI_CALC_TENSOR(
8       I             bi, bj, iMin, iMax, jMin, jMax, K,       I             bi, bj, iMin, iMax, jMin, jMax,
9       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
10       I             myThid )       I             myThid )
11  C     /==========================================================\  C     /==========================================================\
# Line 29  C Line 29  C
29        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
30        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
31        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
32        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER bi,bj,iMin,iMax,jMin,jMax
33        INTEGER myThid        INTEGER myThid
34  CEndOfInterface  CEndOfInterface
35    
36  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
37    
38  C     == Local variables ==  C     == Local variables ==
39        INTEGER i,j,km1,kp1        INTEGER i,j,k,km1,kp1
40        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
42        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
43        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
44        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
45          _RL maskp1, Kgm_tmp
46    
47  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
48        _RS deltaH,zero_rs        _RS deltaH,zero_rs
# Line 50  C     == Local variables == Line 51  C     == Local variables ==
51        _RL Ssq        _RL Ssq
52  #endif  #endif
53    
54          DO k=2,Nr
55        km1=max(1,K-1)  C-- 1rst loop on k : compute Tensor Coeff. at W points.
56        kp1=min(Nr,K)  c      km1 = MAX(1,k-1)
   
57    
58  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
59  !HPF$ INDEPENDENT  !HPF$ INDEPENDENT
# Line 65  C     == Local variables == Line 65  C     == Local variables ==
65         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
66    
67  C      Gradient of Sigma at rVel points  C      Gradient of Sigma at rVel points
68          SlopeX(i,j)=0.25*( sigmaX(i+1, j ,km1) +sigmaX(i,j,km1)          SlopeX(i,j)=0.25*( sigmaX(i+1, j ,k-1) +sigmaX(i,j,k-1)
69       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )
70       &                  *maskC(i,j,k,bi,bj)       &                  *maskC(i,j,k,bi,bj)
71          SlopeY(i,j)=0.25*( sigmaY( i ,j+1,km1) +sigmaY(i,j,km1)          SlopeY(i,j)=0.25*( sigmaY( i ,j+1,k-1) +sigmaY(i,j,k-1)
72       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )
73       &                  *maskC(i,j,k,bi,bj)       &                  *maskC(i,j,k,bi,bj)
74          dSigmaDrReal(i,j)=sigmaR(i,j,k)          dSigmaDrReal(i,j)=sigmaR(i,j,k)
# Line 78  C      Gradient of Sigma at rVel points Line 78  C      Gradient of Sigma at rVel points
78    
79  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
80        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
81       I             dSigmadRReal,       U             dSigmadRReal,
82       I             rF(K),       I             rF(K),
83       U             SlopeX, SlopeY,       U             SlopeX, SlopeY,
84       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
# Line 93  C       Mask Iso-neutral slopes Line 93  C       Mask Iso-neutral slopes
93          SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)          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)  c       Ssq=SlopeX(i,j)*SlopeX(i,j)+SlopeY(i,j)*SlopeY(i,j)
95    
96  C       Components of Redi/GM tensor  C       Components of Redi/GM tensor
97          Kwx(i,j,k,bi,bj)=2.*SlopeX(i,j)*taperFct(i,j)          Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
98          Kwy(i,j,k,bi,bj)=2.*SlopeY(i,j)*taperFct(i,j)          Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
99          Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)          Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
100    
101  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
# Line 125  C       Now we convert deltaH to a non-d Line 125  C       Now we convert deltaH to a non-d
125       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
126          ENDIF          ENDIF
127    
 C       Limit range that KapGM can take  
         VisbeckK(i,j,bi,bj)=  
      &     min(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)  
   
128  #endif /* GM_VISBECK_VARIABLE_K */  #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    
 #ifdef ALLOW_TIMEAVE  
 C--     Time-average  
         GM_Kwx_T(i,j,k,bi,bj)=GM_Kwx_T(i,j,k,bi,bj)  
      &                       +Kwx(i,j,k,bi,bj)*deltaTclock  
         GM_Kwy_T(i,j,k,bi,bj)=GM_Kwy_T(i,j,k,bi,bj)  
      &                       +Kwy(i,j,k,bi,bj)*deltaTclock  
         GM_Kwz_T(i,j,k,bi,bj)=GM_Kwz_T(i,j,k,bi,bj)  
      &                       +Kwz(i,j,k,bi,bj)*deltaTclock  
137  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
138          IF (K.EQ.Nr)        IF ( GM_Visbeck_alpha.NE.0. ) THEN
139       &  Visbeck_K_T(i,j,bi,bj)=Visbeck_K_T(i,j,bi,bj)  C-     Limit range that KapGM can take
140       &                       +VisbeckK(i,j,bi,bj)*deltaTclock         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  #endif
148  #endif /* ALLOW_TIMEAVE */          ENDDO
149         ENDDO         ENDDO
150        ENDDO        ENDIF
151    #endif /* GM_VISBECK_VARIABLE_K */
152    
153  #ifdef ALLOW_TIMEAVE  
154        GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock  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  #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    #ifdef GM_VISBECK_VARIABLE_K
172         &                    + VisbeckK(i,j,bi,bj)
173    #endif
174         &                    )*Kwz(i,j,k,bi,bj)
175           ENDDO
176          ENDDO
177    
178    #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
179    
 #ifdef GM_NON_UNITY_DIAGONAL  
180  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
181        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
182         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
183          SlopeX(i,j)=sigmaX(i,j,km1)          SlopeX(i,j)=sigmaX(i,j,k)
184       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
185          SlopeY(i,j)=0.25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)          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) )       &                    +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )
187       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
188          dSigmaDrReal(i,j)=0.25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )          dSigmaDrReal(i,j)=0.25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )
189       &                          +sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1) )       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )
190       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
191         ENDDO         ENDDO
192        ENDDO        ENDDO
193    
194  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
195        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
196       I             dSigmadRReal,       U             dSigmadRReal,
197       I             rF(K),       I             rF(K),
198       U             SlopeX, SlopeY,       U             SlopeX, SlopeY,
199       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
200       I             bi, bj, myThid )       I             bi, bj, myThid )
201    
202        DO j=1-Oly+1,sNy+Oly-1  #ifdef GM_NON_UNITY_DIAGONAL
203         DO i=1-Olx+1,sNx+Olx-1          DO j=1-Oly+1,sNy+Oly-1
204          Kux(i,j,k,bi,bj)=taperFct(i,j)           DO i=1-Olx+1,sNx+Olx-1
205         ENDDO            Kux(i,j,k,bi,bj) =
206        ENDDO       &     ( 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    
231  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
232        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
# Line 189  C     Gradient of Sigma at V points Line 234  C     Gradient of Sigma at V points
234          SlopeX(i,j)=0.25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)          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) )       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )
236       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
237          SlopeY(i,j)=sigmaY(i,j,km1)          SlopeY(i,j)=sigmaY(i,j,k)
238       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
239          dSigmaDrReal(i,j)=0.25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )          dSigmaDrReal(i,j)=0.25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )
240       &                          +sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1) )       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )
241       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
242         ENDDO         ENDDO
243        ENDDO        ENDDO
244    
245  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
246        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
247       I             dSigmadRReal,       U             dSigmadRReal,
248       I             rF(K),       I             rF(K),
249       U             SlopeX, SlopeY,       U             SlopeX, SlopeY,
250       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
251       I             bi, bj, myThid )       I             bi, bj, myThid )
252    
253    #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        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
287         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
288          Kvy(i,j,k,bi,bj)=taperFct(i,j)          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         ENDDO         ENDDO
295        ENDDO        ENDDO
296          GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
297    #endif /* ALLOW_TIMEAVE */
298    
299  #endif /* GM_NON_UNITY_DIAGONAL */  C-- end 2nd loop on vertical level index k
300          ENDDO
301    
302    
303    #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    
312  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
313    
# Line 222  C     Calculate slopes for use in tensor Line 316  C     Calculate slopes for use in tensor
316    
317    
318        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
319       I             bi, bj, iMin, iMax, jMin, jMax, K,       I             bi, bj, iMin, iMax, jMin, jMax,
320       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
321       I             myThid )       I             myThid )
322  C     /==========================================================\  C     /==========================================================\
# Line 245  C Line 339  C
339        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
340        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
341        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
342        INTEGER bi,bj,iMin,iMax,jMin,jMax,K        INTEGER bi,bj,iMin,iMax,jMin,jMax
343        INTEGER myThid        INTEGER myThid
344  CEndOfInterface  CEndOfInterface
345    
346        INTEGER i, j        INTEGER i, j, k
347    
348  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
349    
350        DO j=1-Oly+1,sNy+Oly-1        DO k=1,Nr
351         DO i=1-Olx+1,sNx+Olx-1         DO j=1-Oly+1,sNy+Oly-1
352          Kwx(i,j,k,bi,bj) = 0.0          DO i=1-Olx+1,sNx+Olx-1
353          Kwy(i,j,k,bi,bj) = 0.0           Kwx(i,j,k,bi,bj) = 0.0
354          Kwz(i,j,k,bi,bj) = 0.0           Kwy(i,j,k,bi,bj) = 0.0
355             Kwz(i,j,k,bi,bj) = 0.0
356            ENDDO
357         ENDDO         ENDDO
358        ENDDO        ENDDO
359  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
360    
361        end        RETURN
362          END
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.22