/[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.11 by mlosch, Wed Sep 25 19:36:50 2002 UTC revision 1.13 by heimbach, Thu Nov 28 17:30:34 2002 UTC
# Line 44  C     == Local variables == Line 44  C     == Local variables ==
44        INTEGER i,j,k,km1,kp1        INTEGER i,j,k,km1,kp1
45        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
46        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47          _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48          _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
# Line 64  C     == Local variables == Line 66  C     == Local variables ==
66            act3 = myThid - 1            act3 = myThid - 1
67            max3 = nTx*nTy            max3 = nTx*nTy
68            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
69            ikey = (act1 + 1) + act2*max1            igmkey = (act1 + 1) + act2*max1
70       &                      + act3*max1*max2       &                      + act3*max1*max2
71       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
72  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
73    
74    #ifdef GM_VISBECK_VARIABLE_K
75          DO j=1-Oly,sNy+Oly
76           DO i=1-Olx,sNx+Olx
77            VisbeckK(i,j,bi,bj) = 0. _d 0
78           ENDDO
79          ENDDO
80    #endif
81    
82        DO k=2,Nr        DO k=2,Nr
83  C-- 1rst loop on k : compute Tensor Coeff. at W points.  C-- 1rst loop on k : compute Tensor Coeff. at W points.
84         km1 = MAX(1,k-1)         km1 = MAX(1,k-1)
# Line 76  C-- 1rst loop on k : compute Tensor Coef Line 86  C-- 1rst loop on k : compute Tensor Coef
86         IF (k.LE.1) maskm1 = 0. _d 0         IF (k.LE.1) maskm1 = 0. _d 0
87    
88  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
89         kkey = (ikey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
90         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
91          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
92           SlopeX(i,j)       = 0. _d 0           SlopeX(i,j)       = 0. _d 0
93           SlopeY(i,j)       = 0. _d 0           SlopeY(i,j)       = 0. _d 0
94             dSigmaDx(i,j)     = 0. _d 0
95             dSigmaDy(i,j)     = 0. _d 0
96           dSigmaDrReal(i,j) = 0. _d 0           dSigmaDrReal(i,j) = 0. _d 0
97           SlopeSqr(i,j)     = 0. _d 0           SlopeSqr(i,j)     = 0. _d 0
98           taperFct(i,j)     = 0. _d 0           taperFct(i,j)     = 0. _d 0
99           Kwx(i,j,k,bi,bj)  = 0. _d 0           Kwx(i,j,k,bi,bj)  = 0. _d 0
100           Kwy(i,j,k,bi,bj)  = 0. _d 0           Kwy(i,j,k,bi,bj)  = 0. _d 0
101           Kwz(i,j,k,bi,bj)  = 0. _d 0           Kwz(i,j,k,bi,bj)  = 0. _d 0
102    # ifdef GM_NON_UNITY_DIAGONAL
103             Kux(i,j,k,bi,bj)  = 0. _d 0
104             Kvy(i,j,k,bi,bj)  = 0. _d 0
105    # endif
106    # ifdef GM_EXTRA_DIAGONAL
107             Kuz(i,j,k,bi,bj)  = 0. _d 0
108             Kvz(i,j,k,bi,bj)  = 0. _d 0
109    # endif
110    # ifdef GM_BOLUS_ADVEC
111             GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
112             GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
113    # endif
114          ENDDO          ENDDO
115         ENDDO         ENDDO
116  #endif  #endif
117    
118        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
119         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
   
120  C      Gradient of Sigma at rVel points  C      Gradient of Sigma at rVel points
121          SlopeX(i,j)=0.25*( sigmaX(i+1, j ,km1) +sigmaX(i,j,km1)          dSigmaDx(i,j)=0.25*( sigmaX(i+1, j ,km1) +sigmaX(i,j,km1)
122       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )
123       &                  *maskC(i,j,k,bi,bj)*maskm1       &                  *maskC(i,j,k,bi,bj)*maskm1
124          SlopeY(i,j)=0.25*( sigmaY( i ,j+1,km1) +sigmaY(i,j,km1)          dSigmaDy(i,j)=0.25*( sigmaY( i ,j+1,km1) +sigmaY(i,j,km1)
125       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )
126       &                  *maskC(i,j,k,bi,bj)*maskm1       &                  *maskC(i,j,k,bi,bj)*maskm1
127          dSigmaDrReal(i,j)=sigmaR(i,j,k)*maskm1          dSigmaDrReal(i,j)=sigmaR(i,j,k)*maskm1
   
128         ENDDO         ENDDO
129        ENDDO        ENDDO
130    
131  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
132  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
133  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
134  CADJ STORE dsigmadrreal(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dsigmadrreal(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
135  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
136    
137  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
138        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
139       U             dSigmadRReal,       U             dSigmadRReal,
140       I             rF(K),       I             rF(K),K,
141       U             SlopeX, SlopeY,       U             SlopeX, SlopeY,
142         U             dSigmaDx, dSigmaDy,
143       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
144       I             bi, bj, myThid )       I             bi, bj, myThid )
145    
# Line 132  C       Mask Iso-neutral slopes Line 155  C       Mask Iso-neutral slopes
155        ENDDO        ENDDO
156    
157  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
158  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)          = comlev1_bibj_k, key=kkey, byte=isbyte
159  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeY(:,:)          = comlev1_bibj_k, key=kkey, byte=isbyte
160  CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)        = comlev1_bibj_k, key=kkey, byte=isbyte
161    CADJ STORE taperFct(:,:)        = comlev1_bibj_k, key=kkey, byte=isbyte
162    #ifdef GM_VISBECK_VARIABLE_K
163    CADJ STORE dSigmaDrReal(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
164    #endif
165  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
166    
167        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
# Line 165  C       Now we convert deltaH to a non-d Line 192  C       Now we convert deltaH to a non-d
192          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
193    
194          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
195          IF (Ssq(i,j).NE.0.) THEN          IF ( Ssq(i,j).NE.0. .AND. dSigmaDrReal(i,j).NE.0. ) THEN
196           N2= -Gravity*recip_RhoConst*dSigmaDrReal(i,j)           N2= -Gravity*recip_RhoConst*dSigmaDrReal(i,j)
197           SN=sqrt(Ssq(i,j)*N2)           SN=sqrt(Ssq(i,j)*N2)
198           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
# Line 182  C-- end 1rst loop on vertical level inde Line 209  C-- end 1rst loop on vertical level inde
209    
210    
211  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
212    #ifdef ALLOW_AUTODIFF_TAMC
213    CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
214    #endif
215        IF ( GM_Visbeck_alpha.NE.0. ) THEN        IF ( GM_Visbeck_alpha.NE.0. ) THEN
216  C-     Limit range that KapGM can take  C-     Limit range that KapGM can take
217         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
# Line 206  C-- 2nd loop on k : compute Tensor Coeff Line 236  C-- 2nd loop on k : compute Tensor Coeff
236         maskp1 = 1. _d 0         maskp1 = 1. _d 0
237         IF (k.GE.Nr) maskp1 = 0. _d 0         IF (k.GE.Nr) maskp1 = 0. _d 0
238    
239    #ifdef ALLOW_AUTODIFF_TAMC
240           kkey = (igmkey-1)*Nr + k
241    #ifdef GM_VISBECK_VARIABLE_K
242    CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj, key=kkey, byte=isbyte
243    CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj, key=kkey, byte=isbyte
244    CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj, key=kkey, byte=isbyte
245    #endif
246    #endif
247    
248  C-    express the Tensor in term of Diffusivity (= m**2 / s )  C-    express the Tensor in term of Diffusivity (= m**2 / s )
249        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
250         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
# Line 222  C-    express the Tensor in term of Diff Line 261  C-    express the Tensor in term of Diff
261       &                    )*Kwz(i,j,k,bi,bj)       &                    )*Kwz(i,j,k,bi,bj)
262         ENDDO         ENDDO
263        ENDDO        ENDDO
264    #ifdef ALLOW_AUTODIFF_TAMC
265    #ifdef GM_VISBECK_VARIABLE_K
266    CADJ STORE VisbeckK(:,:,bi,bj) =
267    CADJ &     comlev1_bibj, key=kkey, byte=isbyte
268    #endif
269    #endif
270    
271  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
272    
273  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
274        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
275         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
276          SlopeX(i,j)=sigmaX(i,j,k)          dSigmaDx(i,j)=sigmaX(i,j,k)
277       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
278          SlopeY(i,j)=0.25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)          dSigmaDy(i,j)=0.25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)
279       &                    +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )       &                      +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )
280       &          *_maskW(i,j,k,bi,bj)       &          *_maskW(i,j,k,bi,bj)
281          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 )
282       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )
# Line 239  C     Gradient of Sigma at U points Line 284  C     Gradient of Sigma at U points
284         ENDDO         ENDDO
285        ENDDO        ENDDO
286    
287    #ifdef ALLOW_AUTODIFF_TAMC
288    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
289    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
290    CADJ STORE dsigmadrreal(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
291    #endif /* ALLOW_AUTODIFF_TAMC */
292    
293  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
294        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
295       U             dSigmadRReal,       U             dSigmadRReal,
296       I             rF(K),       I             rF(K),K,
297       U             SlopeX, SlopeY,       U             SlopeX, SlopeY,
298         U             dSigmaDx, dSigmaDy,
299       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
300       I             bi, bj, myThid )       I             bi, bj, myThid )
301    
# Line 259  C     Calculate slopes for use in tensor Line 311  C     Calculate slopes for use in tensor
311       &     *taperFct(i,j)       &     *taperFct(i,j)
312           ENDDO           ENDDO
313          ENDDO          ENDDO
314    #ifdef ALLOW_AUTODIFF_TAMC
315    # ifndef GM_TAPER_ORIG_CLIPPING
316    CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
317    # endif
318    #endif
319          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
320           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
321            Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )            Kux(i,j,k,bi,bj) = MAX( Kux(i,j,k,bi,bj), GM_Kmin_horiz )
# Line 267  C     Calculate slopes for use in tensor Line 324  C     Calculate slopes for use in tensor
324  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
325    
326  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
327    
328    #ifdef ALLOW_AUTODIFF_TAMC
329    CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
330    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
331    #endif /* ALLOW_AUTODIFF_TAMC */
332        IF (GM_ExtraDiag) THEN        IF (GM_ExtraDiag) THEN
333          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
334           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
# Line 284  C     Calculate slopes for use in tensor Line 346  C     Calculate slopes for use in tensor
346  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
347        DO j=1-Oly+1,sNy+Oly-1        DO j=1-Oly+1,sNy+Oly-1
348         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
349          SlopeX(i,j)=0.25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)          dSigmaDx(i,j)=0.25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
350       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )
351       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
352          SlopeY(i,j)=sigmaY(i,j,k)          dSigmaDy(i,j)=sigmaY(i,j,k)
353       &          *_maskS(i,j,k,bi,bj)       &          *_maskS(i,j,k,bi,bj)
354          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 )
355       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )
# Line 295  C     Gradient of Sigma at V points Line 357  C     Gradient of Sigma at V points
357         ENDDO         ENDDO
358        ENDDO        ENDDO
359    
360    #ifdef ALLOW_AUTODIFF_TAMC
361    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
362    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
363    CADJ STORE dsigmadrreal(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
364    #endif /* ALLOW_AUTODIFF_TAMC */
365    
366  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
367        CALL GMREDI_SLOPE_LIMIT(        CALL GMREDI_SLOPE_LIMIT(
368       U             dSigmadRReal,       U             dSigmadRReal,
369       I             rF(K),       I             rF(K),K,
370       U             SlopeX, SlopeY,       U             SlopeX, SlopeY,
371         U             dSigmaDx, dSigmaDy,
372       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
373       I             bi, bj, myThid )       I             bi, bj, myThid )
374    
# Line 315  C     Calculate slopes for use in tensor Line 384  C     Calculate slopes for use in tensor
384       &     *taperFct(i,j)       &     *taperFct(i,j)
385           ENDDO           ENDDO
386          ENDDO          ENDDO
387    #ifdef ALLOW_AUTODIFF_TAMC
388    # ifndef GM_TAPER_ORIG_CLIPPING
389    CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
390    # endif
391    #endif
392          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
393           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
394            Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )            Kvy(i,j,k,bi,bj) = MAX( Kvy(i,j,k,bi,bj), GM_Kmin_horiz )
# Line 323  C     Calculate slopes for use in tensor Line 397  C     Calculate slopes for use in tensor
397  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
398    
399  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
400    
401    #ifdef ALLOW_AUTODIFF_TAMC
402    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
403    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
404    #endif /* ALLOW_AUTODIFF_TAMC */
405        IF (GM_ExtraDiag) THEN        IF (GM_ExtraDiag) THEN
406          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
407           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22