/[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.26 by jmc, Thu May 24 22:34:38 2007 UTC revision 1.32 by heimbach, Fri Mar 28 18:48:05 2008 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "GMREDI_OPTIONS.h"  #include "GMREDI_OPTIONS.h"
5    #ifdef ALLOW_KPP
6    # include "KPP_OPTIONS.h"
7    #endif
8  #undef OLD_VISBECK_CALC  #undef OLD_VISBECK_CALC
9    
10  CStartOfInterface  CBOP
11    C     !ROUTINE: GMREDI_CALC_TENSOR
12    C     !INTERFACE:
13        SUBROUTINE GMREDI_CALC_TENSOR(        SUBROUTINE GMREDI_CALC_TENSOR(
14       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
15       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
16       I             myThid )       I             bi, bj, myTime, myIter, myThid )
17  C     /==========================================================\  
18  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     !DESCRIPTION: \bv
19  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     *==========================================================*
20  C     |==========================================================|  C     | SUBROUTINE GMREDI_CALC_TENSOR
21  C     \==========================================================/  C     | o Calculate tensor elements for GM/Redi tensor.
22    C     *==========================================================*
23    C     *==========================================================*
24    C     \ev
25    
26    C     !USES:
27        IMPLICIT NONE        IMPLICIT NONE
28    
29  C     == Global variables ==  C     == Global variables ==
# Line 24  C     == Global variables == Line 34  C     == Global variables ==
34  #include "PARAMS.h"  #include "PARAMS.h"
35  #include "GMREDI.h"  #include "GMREDI.h"
36  #include "GMREDI_TAVE.h"  #include "GMREDI_TAVE.h"
37    #ifdef ALLOW_KPP
38    # include "KPP.h"
39    #endif
40    
41  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
42  #include "tamc.h"  #include "tamc.h"
43  #include "tamc_keys.h"  #include "tamc_keys.h"
44  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
45    
46    C     !INPUT/OUTPUT PARAMETERS:
47  C     == Routine arguments ==  C     == Routine arguments ==
48    C     bi, bj    :: tile indices
49    C     myTime    :: Current time in simulation
50    C     myIter    :: Current iteration number in simulation
51    C     myThid    :: My Thread Id. number
52  C  C
53          INTEGER iMin,iMax,jMin,jMax
54        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
55        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
56        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
57        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi, bj
58          _RL     myTime
59          INTEGER myIter
60        INTEGER myThid        INTEGER myThid
61  CEndOfInterface  CEOP
62    
63  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
64    
65    C     !LOCAL VARIABLES:
66  C     == Local variables ==  C     == Local variables ==
67        INTEGER i,j,k,kp1        INTEGER i,j,k,kp1
68        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
# Line 56  C     == Local variables == Line 78  C     == Local variables ==
78        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
79        _RL Cspd, LrhoInf, LrhoSup, fCoriLoc        _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
80    
81          INTEGER kLow_W (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82          INTEGER kLow_S (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
83          _RL locMixLayer(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
84          _RL baseSlope  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
85          _RL hTransLay  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
86          _RL recipLambda(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
87    
88  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
89  #ifdef OLD_VISBECK_CALC  #ifdef OLD_VISBECK_CALC
90        _RL deltaH,zero_rs        _RL deltaH,zero_rs
# Line 111  C---+----1----+----2----+----3----+----4 Line 140  C---+----1----+----2----+----3----+----4
140  #endif  #endif
141    
142  C--   set ldd97_Lrho (for tapering scheme ldd97):  C--   set ldd97_Lrho (for tapering scheme ldd97):
143        IF (GM_taper_scheme.EQ.'ldd97') THEN        IF ( GM_taper_scheme.EQ.'ldd97' .OR.
144         &     GM_taper_scheme.EQ.'fm07' ) THEN
145         Cspd = 2. _d 0         Cspd = 2. _d 0
146         LrhoInf = 15. _d 3         LrhoInf = 15. _d 3
147         LrhoSup = 100. _d 3         LrhoSup = 100. _d 3
# Line 128  C-     Tracer point location (center): Line 158  C-     Tracer point location (center):
158         ENDDO         ENDDO
159  C-     U point location (West):  C-     U point location (West):
160         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
161            kLow_W(1-Olx,j) = 0
162          ldd97_LrhoW(1-Olx,j) = LrhoSup          ldd97_LrhoW(1-Olx,j) = LrhoSup
163          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
164             kLow_W(i,j) = MIN(kLowC(i-1,j,bi,bj),kLowC(i,j,bi,bj))
165           fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))           fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
166           IF (fCoriLoc.NE.0.) THEN           IF (fCoriLoc.NE.0.) THEN
167             ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)             ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
# Line 141  C-     U point location (West): Line 173  C-     U point location (West):
173         ENDDO         ENDDO
174  C-     V point location (South):  C-     V point location (South):
175         DO i=1-Olx+1,sNx+Olx         DO i=1-Olx+1,sNx+Olx
176             kLow_S(i,1-Oly) = 0
177           ldd97_LrhoS(i,1-Oly) = LrhoSup           ldd97_LrhoS(i,1-Oly) = LrhoSup
178         ENDDO         ENDDO
179         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
180          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
181             kLow_S(i,j) = MIN(kLowC(i,j-1,bi,bj),kLowC(i,j,bi,bj))
182           fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))           fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
183           IF (fCoriLoc.NE.0.) THEN           IF (fCoriLoc.NE.0.) THEN
184             ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)             ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
# Line 164  C-    Just initialize to zero (not use a Line 198  C-    Just initialize to zero (not use a
198          ENDDO          ENDDO
199         ENDDO         ENDDO
200        ENDIF        ENDIF
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
201    
202        DO k=2,Nr  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
203  C-- 1rst loop on k : compute Tensor Coeff. at W points.  C-- 1rst loop on k : compute Tensor Coeff. at W points.
204    
205          DO j=1-Oly,sNy+Oly
206           DO i=1-Olx,sNx+Olx
207             hTransLay(i,j) = R_low(i,j,bi,bj)
208             baseSlope(i,j) =  0. _d 0
209             recipLambda(i,j) = 0. _d 0
210             locMixLayer(i,j) = 0. _d 0
211           ENDDO
212          ENDDO
213    #ifdef ALLOW_KPP
214          IF ( useKPP ) THEN
215           DO j=1-Oly,sNy+Oly
216            DO i=1-Olx,sNx+Olx
217             locMixLayer(i,j) = KPPhbl(i,j,bi,bj)
218            ENDDO
219           ENDDO
220          ELSE
221    #else
222          IF ( .TRUE. ) THEN
223    #endif
224           DO j=1-Oly,sNy+Oly
225            DO i=1-Olx,sNx+Olx
226             locMixLayer(i,j) = hMixLayer(i,j,bi,bj)
227            ENDDO
228           ENDDO
229          ENDIF
230    
231          DO k=Nr,2,-1
232    
233  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
234         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
235         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 216  C      Gradient of Sigma at rVel points Line 277  C      Gradient of Sigma at rVel points
277  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
278  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
279  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
280    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
281    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
282    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
283  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
284    
285  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
# Line 276  C     Calculate slopes for use in tensor Line 340  C     Calculate slopes for use in tensor
340         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
341       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
342       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
343         U             hTransLay, baseSlope, recipLambda,
344       U             dSigmaDr,       U             dSigmaDr,
345       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
346       I             ldd97_LrhoC,rF(k),k,       I             ldd97_LrhoC, locMixLayer, rF,
347       I             bi, bj, myThid )       I             kLowC(1-Olx,1-Oly,bi,bj),
348         I             k, bi, bj, myTime, myIter, myThid )
349    
350         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
351          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
# Line 298  CADJ STORE dSigmaDr(:,:)     = comlev1_b Line 364  CADJ STORE dSigmaDr(:,:)     = comlev1_b
364  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
365  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
366    
367    C      Components of Redi/GM tensor
368         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
369          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
370  C      Components of Redi/GM tensor            Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
371           Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)            Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
372           Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)            Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
373           Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)          ENDDO
374           ENDDO
375    
376  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
377  #ifdef OLD_VISBECK_CALC  #ifdef OLD_VISBECK_CALC
378           DO j=1-Oly+1,sNy+Oly-1
379            DO i=1-Olx+1,sNx+Olx-1
380    
381  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
382  C           but do not know if *taperFct (or **2 ?) is necessary  C           but do not know if *taperFct (or **2 ?) is necessary
383          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
384    
385  C--     Depth average of M^2/N^2 * N  C--     Depth average of M^2/N^2 * N
# Line 333  C       Now we convert deltaH to a non-d Line 403  C       Now we convert deltaH to a non-d
403       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
404          ENDIF          ENDIF
405    
 #endif /* OLD_VISBECK_CALC */  
 #endif /* GM_VISBECK_VARIABLE_K */  
406          ENDDO          ENDDO
407         ENDDO         ENDDO
408    #endif /* OLD_VISBECK_CALC */
409    #endif /* GM_VISBECK_VARIABLE_K */
410    
411  C-- end 1rst loop on vertical level index k  C-- end 1rst loop on vertical level index k
412        ENDDO        ENDDO
# Line 346  C-- end 1rst loop on vertical level inde Line 416  C-- end 1rst loop on vertical level inde
416  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
417  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
418  #endif  #endif
419        IF ( GM_Visbeck_alpha.NE.0. ) THEN        IF ( GM_Visbeck_alpha.GT.0. ) THEN
420  C-     Limit range that KapGM can take  C-     Limit range that KapGM can take
421         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
422          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
# Line 362  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1 Line 432  CADJ STORE VisbeckK(:,:,bi,bj) = comlev1
432  cph)  cph)
433  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
434    
435    C-    express the Tensor in term of Diffusivity (= m**2 / s )
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  
   
 C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.  
436        DO k=1,Nr        DO k=1,Nr
        kp1 = MIN(Nr,k+1)  
        maskp1 = 1. _d 0  
        IF (k.GE.Nr) maskp1 = 0. _d 0  
   
437  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
438         kkey = (igmkey-1)*Nr + k         kkey = (igmkey-1)*Nr + k
439  #if (defined (GM_NON_UNITY_DIAGONAL) || \  # if (defined (GM_NON_UNITY_DIAGONAL) || \
440       defined (GM_VISBECK_VARIABLE_K))       defined (GM_VISBECK_VARIABLE_K))
441  CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
442  CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
443  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
444    # endif
445  #endif  #endif
 #endif  
   
 C-    express the Tensor in term of Diffusivity (= m**2 / s )  
446         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
447          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
448  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPREDI_CONTROL
449           Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)           Kgm_tmp = kapredi(i,j,k,bi,bj)
450  #else  #else
451           Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K           Kgm_tmp = GM_isopycK
452    #endif
453    #if (defined (ALLOW_AUTODIFF) && defined (ALLOW_KAPGM_CONTROL))
454         &           + GM_skewflx*kapgm(i,j,k,bi,bj)
455    #else
456         &           + GM_skewflx*GM_background_K
457  #endif  #endif
458  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
459       &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)       &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
460  #endif  #endif
461           Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)           Kwx(i,j,k,bi,bj)= Kgm_tmp*Kwx(i,j,k,bi,bj)
462           Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)           Kwy(i,j,k,bi,bj)= Kgm_tmp*Kwy(i,j,k,bi,bj)
463    #ifdef ALLOW_KAPREDI_CONTROL
464             Kwz(i,j,k,bi,bj)= ( kapredi(i,j,k,bi,bj)
465    #else
466           Kwz(i,j,k,bi,bj)= ( GM_isopycK           Kwz(i,j,k,bi,bj)= ( GM_isopycK
467    #endif
468  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
469       &                     + VisbeckK(i,j,bi,bj)       &                     + VisbeckK(i,j,bi,bj)
470  #endif  #endif
471       &                     )*Kwz(i,j,k,bi,bj)       &                     )*Kwz(i,j,k,bi,bj)
472          ENDDO          ENDDO
473         ENDDO         ENDDO
474          ENDDO
475    
476    #ifdef ALLOW_DIAGNOSTICS
477          IF ( useDiagnostics .AND. GM_taper_scheme.EQ.'fm07' ) THEN
478           CALL DIAGNOSTICS_FILL( hTransLay, 'GM_hTrsL', 0,1,2,bi,bj,myThid)
479           CALL DIAGNOSTICS_FILL( baseSlope, 'GM_baseS', 0,1,2,bi,bj,myThid)
480           CALL DIAGNOSTICS_FILL(recipLambda,'GM_rLamb', 0,1,2,bi,bj,myThid)
481          ENDIF
482    #endif /* ALLOW_DIAGNOSTICS */
483    
484    
485  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
486    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
487    C-- 2nd  k loop : compute Tensor Coeff. at U point
488    
489    #ifdef ALLOW_KPP
490          IF ( useKPP ) THEN
491           DO j=1-Oly,sNy+Oly
492            DO i=2-Olx,sNx+Olx
493             locMixLayer(i,j) = ( KPPhbl(i-1,j,bi,bj)
494         &                      + KPPhbl( i ,j,bi,bj) )*op5
495            ENDDO
496           ENDDO
497          ELSE
498    #else
499          IF ( .TRUE. ) THEN
500    #endif
501           DO j=1-Oly,sNy+Oly
502            DO i=2-Olx,sNx+Olx
503             locMixLayer(i,j) = ( hMixLayer(i-1,j,bi,bj)
504         &                      + hMixLayer( i ,j,bi,bj) )*op5
505            ENDDO
506           ENDDO
507          ENDIF
508          DO j=1-Oly,sNy+Oly
509           DO i=1-Olx,sNx+Olx
510             hTransLay(i,j) =  0.
511             baseSlope(i,j) =  0.
512             recipLambda(i,j)= 0.
513           ENDDO
514           DO i=2-Olx,sNx+Olx
515             hTransLay(i,j) = MAX( R_low(i-1,j,bi,bj), R_low(i,j,bi,bj) )
516           ENDDO
517          ENDDO
518    
519          DO k=Nr,1,-1
520           kp1 = MIN(Nr,k+1)
521           maskp1 = 1. _d 0
522           IF (k.GE.Nr) maskp1 = 0. _d 0
523    #ifdef ALLOW_AUTODIFF_TAMC
524           kkey = (igmkey-1)*Nr + k
525    #endif
526    
527  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
528         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
# Line 422  C     Gradient of Sigma at U points Line 542  C     Gradient of Sigma at U points
542  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
543  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
544  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
545  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
546    CADJ STORE locMixLayer(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
547    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
548    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
549    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
550  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
551    
552  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
553         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
554       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
555       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
556         U             hTransLay, baseSlope, recipLambda,
557       U             dSigmaDr,       U             dSigmaDr,
558       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
559       I             ldd97_LrhoW,rC(k),k,       I             ldd97_LrhoW, locMixLayer, rC,
560       I             bi, bj, myThid )       I             kLow_W,
561         I             k, bi, bj, myTime, myIter, myThid )
562    
 cph( NEW  
563  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 cph(  
564  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
565  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
 cph)  
566  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 cph)  
567    
568  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
569  c      IF ( GM_nonUnitDiag ) THEN  c      IF ( GM_nonUnitDiag ) THEN
570          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
571           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
572            Kux(i,j,k,bi,bj) =            Kux(i,j,k,bi,bj) =
573    #ifdef ALLOW_KAPREDI_CONTROL
574         &     ( kapredi(i,j,k,bi,bj)
575    #else
576       &     ( GM_isopycK       &     ( GM_isopycK
577    #endif
578  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
579       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
580  #endif  #endif
581       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
582           ENDDO           ENDDO
583          ENDDO          ENDDO
584  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 475  c      ENDIF Line 600  c      ENDIF
600  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
601  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
602  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
603         IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
604          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
605           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
606            Kuz(i,j,k,bi,bj) =            Kuz(i,j,k,bi,bj) =
607  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPREDI_CONTROL
608       &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)       &     ( kapredi(i,j,k,bi,bj)
609  #else  #else
610       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     ( GM_isopycK
611    #endif
612    #if (defined (ALLOW_AUTODIFF) && defined (ALLOW_KAPGM_CONTROL))
613         &     - GM_skewflx*kapgm(i,j,k,bi,bj)
614    #else
615         &     - GM_skewflx*GM_background_K
616  #endif  #endif
617  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
618       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*GM_advect
# Line 499  CADJ STORE taperFct(:,:)     = comlev1_b Line 629  CADJ STORE taperFct(:,:)     = comlev1_b
629          DO j=1,sNy          DO j=1,sNy
630           DO i=1,sNx+1           DO i=1,sNx+1
631  C         store in tmp1k Kuz_Redi  C         store in tmp1k Kuz_Redi
632    #ifdef ALLOW_KAPREDI_CONTROL
633              tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
634    #else
635            tmp1k(i,j) = ( GM_isopycK            tmp1k(i,j) = ( GM_isopycK
636    #endif
637  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
638       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
639  #endif  #endif
# Line 531  C-        Vertical gradients interpolate Line 665  C-        Vertical gradients interpolate
665         ENDIF         ENDIF
666  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
667    
668    C-- end 2nd  loop on vertical level index k
669          ENDDO
670    
671    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
672    C-- 3rd  k loop : compute Tensor Coeff. at V point
673    
674    #ifdef ALLOW_KPP
675          IF ( useKPP ) THEN
676           DO j=2-Oly,sNy+Oly
677            DO i=1-Olx,sNx+Olx
678             locMixLayer(i,j) = ( KPPhbl(i,j-1,bi,bj)
679         &                      + KPPhbl(i, j ,bi,bj) )*op5
680            ENDDO
681           ENDDO
682          ELSE
683    #else
684          IF ( .TRUE. ) THEN
685    #endif
686           DO j=2-Oly,sNy+Oly
687            DO i=1-Olx,sNx+Olx
688             locMixLayer(i,j) = ( hMixLayer(i,j-1,bi,bj)
689         &                      + hMixLayer(i, j ,bi,bj) )*op5
690            ENDDO
691           ENDDO
692          ENDIF
693          DO j=1-Oly,sNy+Oly
694           DO i=1-Olx,sNx+Olx
695             hTransLay(i,j) =  0.
696             baseSlope(i,j) =  0.
697             recipLambda(i,j)= 0.
698           ENDDO
699          ENDDO
700          DO j=2-Oly,sNy+Oly
701           DO i=1-Olx,sNx+Olx
702             hTransLay(i,j) = MAX( R_low(i,j-1,bi,bj), R_low(i,j,bi,bj) )
703           ENDDO
704          ENDDO
705    
706  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
707          DO k=Nr,1,-1
708           kp1 = MIN(Nr,k+1)
709           maskp1 = 1. _d 0
710           IF (k.GE.Nr) maskp1 = 0. _d 0
711    #ifdef ALLOW_AUTODIFF_TAMC
712           kkey = (igmkey-1)*Nr + k
713    #endif
714    
715         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
716          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
717           dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)           dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
# Line 548  C     Gradient of Sigma at V points Line 728  C     Gradient of Sigma at V points
728  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
729  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
730  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
731  CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
732    CADJ STORE baseSlope(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
733    CADJ STORE hTransLay(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
734    CADJ STORE recipLambda(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
735  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
736    
737  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
738         CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
739       O             SlopeX, SlopeY,       O             SlopeX, SlopeY,
740       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
741         U             hTransLay, baseSlope, recipLambda,
742       U             dSigmaDr,       U             dSigmaDr,
743       I             dSigmaDx, dSigmaDy,       I             dSigmaDx, dSigmaDy,
744       I             ldd97_LrhoS,rC(k),k,       I             ldd97_LrhoS, locMixLayer, rC,
745       I             bi, bj, myThid )       I             kLow_S,
746         I             k, bi, bj, myTime, myIter, myThid )
747    
748  cph(  cph(
749  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 573  c      IF ( GM_nonUnitDiag ) THEN Line 758  c      IF ( GM_nonUnitDiag ) THEN
758          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
759           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
760            Kvy(i,j,k,bi,bj) =            Kvy(i,j,k,bi,bj) =
761    #ifdef ALLOW_KAPREDI_CONTROL
762         &     ( kapredi(i,j,k,bi,bj)
763    #else
764       &     ( GM_isopycK       &     ( GM_isopycK
765    #endif
766  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
767       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
768  #endif  #endif
769       &     )       &     )*taperFct(i,j)
      &     *taperFct(i,j)  
770           ENDDO           ENDDO
771          ENDDO          ENDDO
772  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 600  c      ENDIF Line 788  c      ENDIF
788  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
789  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
790  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
791         IF (GM_ExtraDiag) THEN         IF ( GM_ExtraDiag ) THEN
792          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
793           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
794            Kvz(i,j,k,bi,bj) =            Kvz(i,j,k,bi,bj) =
795  #ifdef ALLOW_KAPGM_CONTROL  #ifdef ALLOW_KAPREDI_CONTROL
796       &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)       &     ( kapredi(i,j,k,bi,bj)
797    #else
798         &     ( GM_isopycK
799    #endif
800    #if (defined (ALLOW_AUTODIFF) && defined (ALLOW_KAPGM_CONTROL))
801         &     - GM_skewflx*kapgm(i,j,k,bi,bj)
802  #else  #else
803       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     - GM_skewflx*GM_background_K
804  #endif  #endif
805  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
806       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*GM_advect
# Line 620  CADJ STORE taperFct(:,:)     = comlev1_b Line 813  CADJ STORE taperFct(:,:)     = comlev1_b
813    
814  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
815         IF (doDiagRediFlx) THEN         IF (doDiagRediFlx) THEN
816  c       km1 = MAX(k-1,1)          km1 = MAX(k-1,1)
817          DO j=1,sNy+1          DO j=1,sNy+1
818           DO i=1,sNx           DO i=1,sNx
819  C         store in tmp1k Kvz_Redi  C         store in tmp1k Kvz_Redi
820    #ifdef ALLOW_KAPREDI_CONTROL
821              tmp1k(i,j) = ( kapredi(i,j,k,bi,bj)
822    #else
823            tmp1k(i,j) = ( GM_isopycK            tmp1k(i,j) = ( GM_isopycK
824    #endif
825  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
826       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0       &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
827  #endif  #endif
# Line 656  C-        Vertical gradients interpolate Line 853  C-        Vertical gradients interpolate
853         ENDIF         ENDIF
854  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
855    
856  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */  C-- end 3rd  loop on vertical level index k
   
 C-- end 2nd loop on vertical level index k  
857        ENDDO        ENDDO
858    
859    #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
860    
861    
862  #ifdef GM_BOLUS_ADVEC  #ifdef GM_BOLUS_ADVEC
863        IF (GM_AdvForm) THEN        IF (GM_AdvForm) THEN
# Line 716  C--   Time-average Line 913  C--   Time-average
913    
914    
915        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(        SUBROUTINE GMREDI_CALC_TENSOR_DUMMY(
916       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
917       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
918       I             myThid )       I             bi, bj, myTime, myIter, myThid )
919  C     /==========================================================\  C     /==========================================================\
920  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |  C     | SUBROUTINE GMREDI_CALC_TENSOR                            |
921  C     | o Calculate tensor elements for GM/Redi tensor.          |  C     | o Calculate tensor elements for GM/Redi tensor.          |
# Line 736  C Line 933  C
933        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
934        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
935        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
936        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
937          INTEGER bi, bj
938          _RL     myTime
939          INTEGER myIter
940        INTEGER myThid        INTEGER myThid
941  CEndOfInterface  CEndOfInterface
942    
       INTEGER i, j, k  
   
943  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
944    
945          INTEGER i, j, k
946    
947        DO k=1,Nr        DO k=1,Nr
948         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
949          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
Legend:
Removed from v.1.26  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22