/[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.9 by jmc, Sun Dec 16 18:54:49 2001 UTC revision 1.26 by jmc, Thu May 24 22:34:38 2007 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    #undef OLD_VISBECK_CALC
6    
7  CStartOfInterface  CStartOfInterface
8        SUBROUTINE GMREDI_CALC_TENSOR(        SUBROUTINE GMREDI_CALC_TENSOR(
# Line 22  C     == Global variables == Line 23  C     == Global variables ==
23  #include "EEPARAMS.h"  #include "EEPARAMS.h"
24  #include "PARAMS.h"  #include "PARAMS.h"
25  #include "GMREDI.h"  #include "GMREDI.h"
26  #include "GMREDI_DIAGS.h"  #include "GMREDI_TAVE.h"
27    
28    #ifdef ALLOW_AUTODIFF_TAMC
29    #include "tamc.h"
30    #include "tamc_keys.h"
31    #endif /* ALLOW_AUTODIFF_TAMC */
32    
33  C     == Routine arguments ==  C     == Routine arguments ==
34  C  C
# Line 36  CEndOfInterface Line 42  CEndOfInterface
42  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
43    
44  C     == Local variables ==  C     == Local variables ==
45        INTEGER i,j,k,km1,kp1        INTEGER i,j,k,kp1
46        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
48        _RL dSigmaDrReal(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
49          _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50          _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL SlopeSqr(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
52        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL taperFct(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
53        _RL maskp1, Kgm_tmp        _RL maskp1, Kgm_tmp
54          _RL ldd97_LrhoC(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
55          _RL ldd97_LrhoW(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56          _RL ldd97_LrhoS(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57          _RL Cspd, LrhoInf, LrhoSup, fCoriLoc
58    
59  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
60        _RS deltaH,zero_rs  #ifdef OLD_VISBECK_CALC
61        PARAMETER(zero_rs=0.)        _RL deltaH,zero_rs
62          PARAMETER(zero_rs=0.D0)
63        _RL N2,SN        _RL N2,SN
64        _RL Ssq        _RL Ssq(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
65    #else
66          _RL dSigmaH
67          _RL deltaH, integrDepth
68          _RL Sloc, M2loc, SNloc
69    #endif
70  #endif  #endif
71    
72        DO k=2,Nr  #ifdef ALLOW_DIAGNOSTICS
73  C-- 1rst loop on k : compute Tensor Coeff. at W points.        LOGICAL  doDiagRediFlx
74  c      km1 = MAX(1,k-1)        LOGICAL  DIAGNOSTICS_IS_ON
75          EXTERNAL DIAGNOSTICS_IS_ON
76          INTEGER  km1
77          _RL dTdz
78          _RL tmp1k(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79    #endif
80    
81    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
82    
83  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
84  !HPF$ INDEPENDENT            act1 = bi - myBxLo(myThid)
85              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
86              act2 = bj - myByLo(myThid)
87              max2 = myByHi(myThid) - myByLo(myThid) + 1
88              act3 = myThid - 1
89              max3 = nTx*nTy
90              act4 = ikey_dynamics - 1
91              igmkey = (act1 + 1) + act2*max1
92         &                      + act3*max1*max2
93         &                      + act4*max1*max2*max3
94    #endif /* ALLOW_AUTODIFF_TAMC */
95    
96    #ifdef ALLOW_DIAGNOSTICS
97          doDiagRediFlx = .FALSE.
98          IF ( useDiagnostics ) THEN
99            doDiagRediFlx = DIAGNOSTICS_IS_ON('GM_KuzTz', myThid )
100            doDiagRediFlx = doDiagRediFlx .OR.
101         &                  DIAGNOSTICS_IS_ON('GM_KvzTz', myThid )
102          ENDIF
103  #endif  #endif
104        DO j=1-Oly+1,sNy+Oly-1  
105    #ifdef GM_VISBECK_VARIABLE_K
106          DO j=1-Oly,sNy+Oly
107           DO i=1-Olx,sNx+Olx
108            VisbeckK(i,j,bi,bj) = 0. _d 0
109           ENDDO
110          ENDDO
111    #endif
112    
113    C--   set ldd97_Lrho (for tapering scheme ldd97):
114          IF (GM_taper_scheme.EQ.'ldd97') THEN
115           Cspd = 2. _d 0
116           LrhoInf = 15. _d 3
117           LrhoSup = 100. _d 3
118    C-     Tracer point location (center):
119           DO j=1-Oly,sNy+Oly
120            DO i=1-Olx,sNx+Olx
121             IF (fCori(i,j,bi,bj).NE.0.) THEN
122               ldd97_LrhoC(i,j) = Cspd/ABS(fCori(i,j,bi,bj))
123             ELSE
124               ldd97_LrhoC(i,j) = LrhoSup
125             ENDIF
126             ldd97_LrhoC(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoC(i,j),LrhoSup))
127            ENDDO
128           ENDDO
129    C-     U point location (West):
130           DO j=1-Oly,sNy+Oly
131            ldd97_LrhoW(1-Olx,j) = LrhoSup
132            DO i=1-Olx+1,sNx+Olx
133             fCoriLoc = op5*(fCori(i-1,j,bi,bj)+fCori(i,j,bi,bj))
134             IF (fCoriLoc.NE.0.) THEN
135               ldd97_LrhoW(i,j) = Cspd/ABS(fCoriLoc)
136             ELSE
137               ldd97_LrhoW(i,j) = LrhoSup
138             ENDIF
139             ldd97_LrhoW(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoW(i,j),LrhoSup))
140            ENDDO
141           ENDDO
142    C-     V point location (South):
143           DO i=1-Olx+1,sNx+Olx
144             ldd97_LrhoS(i,1-Oly) = LrhoSup
145           ENDDO
146           DO j=1-Oly+1,sNy+Oly
147            DO i=1-Olx,sNx+Olx
148             fCoriLoc = op5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
149             IF (fCoriLoc.NE.0.) THEN
150               ldd97_LrhoS(i,j) = Cspd/ABS(fCoriLoc)
151             ELSE
152               ldd97_LrhoS(i,j) = LrhoSup
153             ENDIF
154             ldd97_LrhoS(i,j) = MAX(LrhoInf,MIN(ldd97_LrhoS(i,j),LrhoSup))
155            ENDDO
156           ENDDO
157          ELSE
158    C-    Just initialize to zero (not use anyway)
159           DO j=1-Oly,sNy+Oly
160            DO i=1-Olx,sNx+Olx
161              ldd97_LrhoC(i,j) = 0. _d 0
162              ldd97_LrhoW(i,j) = 0. _d 0
163              ldd97_LrhoS(i,j) = 0. _d 0
164            ENDDO
165           ENDDO
166          ENDIF
167    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
168    
169          DO k=2,Nr
170    C-- 1rst loop on k : compute Tensor Coeff. at W points.
171    
172  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
173  !HPF$ INDEPENDENT         kkey = (igmkey-1)*Nr + k
174           DO j=1-Oly,sNy+Oly
175            DO i=1-Olx,sNx+Olx
176             SlopeX(i,j)       = 0. _d 0
177             SlopeY(i,j)       = 0. _d 0
178             dSigmaDx(i,j)     = 0. _d 0
179             dSigmaDy(i,j)     = 0. _d 0
180             dSigmaDr(i,j)     = 0. _d 0
181             SlopeSqr(i,j)     = 0. _d 0
182             taperFct(i,j)     = 0. _d 0
183             Kwx(i,j,k,bi,bj)  = 0. _d 0
184             Kwy(i,j,k,bi,bj)  = 0. _d 0
185             Kwz(i,j,k,bi,bj)  = 0. _d 0
186    # ifdef GM_NON_UNITY_DIAGONAL
187             Kux(i,j,k,bi,bj)  = 0. _d 0
188             Kvy(i,j,k,bi,bj)  = 0. _d 0
189    # endif
190    # ifdef GM_EXTRA_DIAGONAL
191             Kuz(i,j,k,bi,bj)  = 0. _d 0
192             Kvz(i,j,k,bi,bj)  = 0. _d 0
193    # endif
194    # ifdef GM_BOLUS_ADVEC
195             GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
196             GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
197    # endif
198            ENDDO
199           ENDDO
200  #endif  #endif
        DO i=1-Olx+1,sNx+Olx-1  
201    
202           DO j=1-Oly+1,sNy+Oly-1
203            DO i=1-Olx+1,sNx+Olx-1
204  C      Gradient of Sigma at rVel points  C      Gradient of Sigma at rVel points
205          SlopeX(i,j)=0.25*( sigmaX(i+1, j ,k-1) +sigmaX(i,j,k-1)           dSigmaDx(i,j)=op25*( sigmaX(i+1,j,k-1)+sigmaX(i,j,k-1)
206       &                    +sigmaX(i+1, j , k ) +sigmaX(i,j, k ) )       &                       +sigmaX(i+1,j, k )+sigmaX(i,j, k )
207       &                  *maskC(i,j,k,bi,bj)       &                      )*maskC(i,j,k,bi,bj)
208          SlopeY(i,j)=0.25*( sigmaY( i ,j+1,k-1) +sigmaY(i,j,k-1)           dSigmaDy(i,j)=op25*( sigmaY(i,j+1,k-1)+sigmaY(i,j,k-1)
209       &                    +sigmaY( i ,j+1, k ) +sigmaY(i,j, k ) )       &                       +sigmaY(i,j+1, k )+sigmaY(i,j, k )
210       &                  *maskC(i,j,k,bi,bj)       &                      )*maskC(i,j,k,bi,bj)
211          dSigmaDrReal(i,j)=sigmaR(i,j,k)           dSigmaDr(i,j)=sigmaR(i,j,k)
212            ENDDO
213         ENDDO         ENDDO
214        ENDDO  
215    #ifdef ALLOW_AUTODIFF_TAMC
216    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
217    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
218    CADJ STORE dSigmaDr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
219    #endif /* ALLOW_AUTODIFF_TAMC */
220    
221    #ifdef GM_VISBECK_VARIABLE_K
222    #ifndef OLD_VISBECK_CALC
223           IF ( GM_Visbeck_alpha.GT.0. .AND.
224         &      -rC(k-1).LT.GM_Visbeck_depth ) THEN
225    
226    C--     Depth average of f/sqrt(Ri) = M^2/N^2 * N
227    C       M^2 and N^2 are horizontal & vertical gradient of buoyancy.
228    
229    C       Calculate terms for mean Richardson number which is used
230    C       in the "variable K" parameterisaton:
231    C       compute depth average from surface down to the bottom or
232    C       GM_Visbeck_depth, whatever is the shallower.
233    
234            DO j=1-Oly+1,sNy+Oly-1
235             DO i=1-Olx+1,sNx+Olx-1
236              IF ( maskC(i,j,k,bi,bj).NE.0. ) THEN
237               integrDepth = -rC( kLowC(i,j,bi,bj) )
238    C-      in 2 steps to avoid mix of RS & RL type in min fct. arguments
239               integrDepth = MIN( integrDepth, GM_Visbeck_depth )
240    C       Distance between level center above and the integration depth
241               deltaH = integrDepth + rC(k-1)
242    C       If negative then we are below the integration level
243    C       (cannot be the case with 2 conditions on maskC & -rC(k-1))
244    C       If positive we limit this to the distance from center above
245               deltaH = MIN( deltaH, drC(k) )
246    C       Now we convert deltaH to a non-dimensional fraction
247               deltaH = deltaH/( integrDepth+rC(1) )
248    
249    C--      compute: ( M^2 * S )^1/2   (= M^2 / N since S=M^2/N^2 )
250               dSigmaH = dSigmaDx(i,j)*dSigmaDx(i,j)
251         &             + dSigmaDy(i,j)*dSigmaDy(i,j)
252               IF ( dSigmaH .GT. 0. _d 0 ) THEN
253                 dSigmaH = SQRT( dSigmaH )
254    C-       compute slope, limited by GM_maxSlope:
255                 IF ( -dSigmaDr(i,j).GT.dSigmaH*GM_rMaxSlope ) THEN
256                  Sloc = dSigmaH / ( -dSigmaDr(i,j) )
257                 ELSE
258                  Sloc = GM_maxSlope
259                 ENDIF
260                 M2loc = Gravity*recip_RhoConst*dSigmaH
261                 SNloc = SQRT( Sloc*M2loc )
262               ELSE
263                 SNloc = 0. _d 0
264               ENDIF
265               VisbeckK(i,j,bi,bj) = VisbeckK(i,j,bi,bj)
266         &       +deltaH*GM_Visbeck_alpha
267         &              *GM_Visbeck_length*GM_Visbeck_length*SNloc
268              ENDIF
269             ENDDO
270            ENDDO
271           ENDIF
272    #endif /* ndef OLD_VISBECK_CALC */
273    #endif /* GM_VISBECK_VARIABLE_K */
274    
275  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
276        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
277       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),  
      U             SlopeX, SlopeY,  
278       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
279         U             dSigmaDr,
280         I             dSigmaDx, dSigmaDy,
281         I             ldd97_LrhoC,rF(k),k,
282       I             bi, bj, myThid )       I             bi, bj, myThid )
283    
284        DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
285         DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
286    C      Mask Iso-neutral slopes
287             SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)
288             SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)
289             SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)
290            ENDDO
291           ENDDO
292    
293    #ifdef ALLOW_AUTODIFF_TAMC
294    CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
295    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
296    CADJ STORE SlopeSqr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
297    CADJ STORE dSigmaDr(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
298    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
299    #endif /* ALLOW_AUTODIFF_TAMC */
300    
301  C       Mask Iso-neutral slopes         DO j=1-Oly+1,sNy+Oly-1
302          SlopeX(i,j)=SlopeX(i,j)*maskC(i,j,k,bi,bj)          DO i=1-Olx+1,sNx+Olx-1
303          SlopeY(i,j)=SlopeY(i,j)*maskC(i,j,k,bi,bj)  C      Components of Redi/GM tensor
304          SlopeSqr(i,j)=SlopeSqr(i,j)*maskC(i,j,k,bi,bj)           Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)
305  c       Ssq=SlopeX(i,j)*SlopeX(i,j)+SlopeY(i,j)*SlopeY(i,j)           Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)
306             Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)
 C       Components of Redi/GM tensor  
         Kwx(i,j,k,bi,bj)= SlopeX(i,j)*taperFct(i,j)  
         Kwy(i,j,k,bi,bj)= SlopeY(i,j)*taperFct(i,j)  
         Kwz(i,j,k,bi,bj)= SlopeSqr(i,j)*taperFct(i,j)  
307    
308  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
309    #ifdef OLD_VISBECK_CALC
310    
311  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K  C- note (jmc) : moved here since only used in VISBECK_VARIABLE_K
312  C           but don't know if *taperFct (or **2 ?) is necessary  C           but do not know if *taperFct (or **2 ?) is necessary
313          Ssq=SlopeSqr(i,j)*taperFct(i,j)          Ssq(i,j)=SlopeSqr(i,j)*taperFct(i,j)
314    
315  C--     Depth average of M^2/N^2 * N  C--     Depth average of M^2/N^2 * N
316    
# Line 118  C       Now we convert deltaH to a non-d Line 326  C       Now we convert deltaH to a non-d
326          deltaH=deltaH/GM_Visbeck_depth          deltaH=deltaH/GM_Visbeck_depth
327    
328          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.          IF (K.eq.2) VisbeckK(i,j,bi,bj)=0.
329          IF (Ssq.NE.0.) THEN          IF ( Ssq(i,j).NE.0. .AND. dSigmaDr(i,j).NE.0. ) THEN
330           N2= -Gravity*recip_Rhonil*dSigmaDrReal(i,j)           N2= -Gravity*recip_RhoConst*dSigmaDr(i,j)
331           SN=sqrt(Ssq*N2)           SN=sqrt(Ssq(i,j)*N2)
332           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH           VisbeckK(i,j,bi,bj)=VisbeckK(i,j,bi,bj)+deltaH
333       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN       &      *GM_Visbeck_alpha*GM_Visbeck_length*GM_Visbeck_length*SN
334          ENDIF          ENDIF
335    
336    #endif /* OLD_VISBECK_CALC */
337  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
338            ENDDO
339         ENDDO         ENDDO
       ENDDO  
340    
341  C-- end 1rst loop on vertical level index k  C-- end 1rst loop on vertical level index k
342        ENDDO        ENDDO
343    
344    
345  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
346    #ifdef ALLOW_AUTODIFF_TAMC
347    CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
348    #endif
349        IF ( GM_Visbeck_alpha.NE.0. ) THEN        IF ( GM_Visbeck_alpha.NE.0. ) THEN
350  C-     Limit range that KapGM can take  C-     Limit range that KapGM can take
351         DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
352          DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
353           VisbeckK(i,j,bi,bj)=           VisbeckK(i,j,bi,bj)=
354       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)       &       MIN(VisbeckK(i,j,bi,bj),GM_Visbeck_maxval_K)
 #ifdef ALLOW_TIMEAVE  
          Visbeck_K_T(i,j,bi,bj)=Visbeck_K_T(i,j,bi,bj)  
      &                         +VisbeckK(i,j,bi,bj)*deltaTclock  
 #endif  
355          ENDDO          ENDDO
356         ENDDO         ENDDO
357        ENDIF        ENDIF
358    cph( NEW
359    #ifdef ALLOW_AUTODIFF_TAMC
360    CADJ STORE VisbeckK(:,:,bi,bj) = comlev1_bibj, key=igmkey, byte=isbyte
361    #endif
362    cph)
363  #endif /* GM_VISBECK_VARIABLE_K */  #endif /* GM_VISBECK_VARIABLE_K */
364    
365    
366  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
367    
368  C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.  C-- 2nd loop on k : compute Tensor Coeff. at U,V levels.
369        DO k=1,Nr        DO k=1,Nr
370         kp1 = MIN(Nr,k+1)         kp1 = MIN(Nr,k+1)
371         maskp1 = 1. _d 0         maskp1 = 1. _d 0
372         IF (k.GE.Nr) maskp1 = 0. _d 0         IF (k.GE.Nr) maskp1 = 0. _d 0
373    
374    #ifdef ALLOW_AUTODIFF_TAMC
375           kkey = (igmkey-1)*Nr + k
376    #if (defined (GM_NON_UNITY_DIAGONAL) || \
377         defined (GM_VISBECK_VARIABLE_K))
378    CADJ STORE Kwx(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
379    CADJ STORE Kwy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
380    CADJ STORE Kwz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
381    #endif
382    #endif
383    
384  C-    express the Tensor in term of Diffusivity (= m**2 / s )  C-    express the Tensor in term of Diffusivity (= m**2 / s )
385        DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
386         DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
387          Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K  #ifdef ALLOW_KAPGM_CONTROL
388             Kgm_tmp = GM_isopycK + GM_skewflx*kapgm(i,j,k,bi,bj)
389    #else
390             Kgm_tmp = GM_isopycK + GM_skewflx*GM_background_K
391    #endif
392  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
393       &          + VisbeckK(i,j,bi,bj)*(1.+GM_skewflx)           &           + VisbeckK(i,j,bi,bj)*(1. _d 0 + GM_skewflx)
394  #endif  #endif
395          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)
396          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)
397          Kwz(i,j,k,bi,bj)= ( GM_isopycK           Kwz(i,j,k,bi,bj)= ( GM_isopycK
398  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
399       &                    + VisbeckK(i,j,bi,bj)       &                     + VisbeckK(i,j,bi,bj)
400  #endif  #endif
401       &                    )*Kwz(i,j,k,bi,bj)       &                     )*Kwz(i,j,k,bi,bj)
402            ENDDO
403         ENDDO         ENDDO
       ENDDO  
404    
405  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )  #if ( defined (GM_NON_UNITY_DIAGONAL) || defined (GM_EXTRA_DIAGONAL) )
406    
407  C     Gradient of Sigma at U points  C     Gradient of Sigma at U points
408        DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
409         DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
410          SlopeX(i,j)=sigmaX(i,j,k)           dSigmaDx(i,j)=sigmaX(i,j,k)
411       &          *_maskW(i,j,k,bi,bj)       &                       *_maskW(i,j,k,bi,bj)
412          SlopeY(i,j)=0.25*( sigmaY(i-1,j+1,k) +sigmaY(i,j+1,k)           dSigmaDy(i,j)=op25*( sigmaY(i-1,j+1,k)+sigmaY(i,j+1,k)
413       &                    +sigmaY(i-1, j ,k) +sigmaY(i, j ,k) )       &                       +sigmaY(i-1, j ,k)+sigmaY(i, j ,k)
414       &          *_maskW(i,j,k,bi,bj)       &                      )*_maskW(i,j,k,bi,bj)
415          dSigmaDrReal(i,j)=0.25*( sigmaR(i-1,j, k ) +sigmaR(i,j, k )           dSigmaDr(i,j)=op25*( sigmaR(i-1,j, k )+sigmaR(i,j, k )
416       &                  +maskp1*(sigmaR(i-1,j,kp1) +sigmaR(i,j,kp1)) )       &                      +(sigmaR(i-1,j,kp1)+sigmaR(i,j,kp1))*maskp1
417       &          *_maskW(i,j,k,bi,bj)       &                      )*_maskW(i,j,k,bi,bj)
418            ENDDO
419         ENDDO         ENDDO
420        ENDDO  
421    #ifdef ALLOW_AUTODIFF_TAMC
422    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
423    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
424    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
425    CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
426    #endif /* ALLOW_AUTODIFF_TAMC */
427    
428  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
429        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
430       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),  
      U             SlopeX, SlopeY,  
431       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
432         U             dSigmaDr,
433         I             dSigmaDx, dSigmaDy,
434         I             ldd97_LrhoW,rC(k),k,
435       I             bi, bj, myThid )       I             bi, bj, myThid )
436    
437    cph( NEW
438    #ifdef ALLOW_AUTODIFF_TAMC
439    cph(
440    CADJ STORE SlopeSqr(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
441    CADJ STORE taperFct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
442    cph)
443    #endif /* ALLOW_AUTODIFF_TAMC */
444    cph)
445    
446  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
447    c      IF ( GM_nonUnitDiag ) THEN
448          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
449           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
450            Kux(i,j,k,bi,bj) =            Kux(i,j,k,bi,bj) =
451       &     ( GM_isopycK       &     ( GM_isopycK
452  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
453       &     +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))
454  #endif  #endif
455       &     )*taperFct(i,j)       &     )
456         &     *taperFct(i,j)
457             ENDDO
458            ENDDO
459    #ifdef ALLOW_AUTODIFF_TAMC
460    # ifdef GM_EXCLUDE_CLIPPING
461    CADJ STORE Kux(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
462    # endif
463    #endif
464            DO j=1-Oly+1,sNy+Oly-1
465             DO i=1-Olx+1,sNx+Olx-1
466            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 )
467           ENDDO           ENDDO
468          ENDDO          ENDDO
469    c      ENDIF
470  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
471    
472  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
473        IF (GM_ExtraDiag) THEN  
474    #ifdef ALLOW_AUTODIFF_TAMC
475    CADJ STORE SlopeX(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
476    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
477    #endif /* ALLOW_AUTODIFF_TAMC */
478           IF (GM_ExtraDiag) THEN
479          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
480           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
481            Kuz(i,j,k,bi,bj) =            Kuz(i,j,k,bi,bj) =
482    #ifdef ALLOW_KAPGM_CONTROL
483         &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
484    #else
485       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     ( GM_isopycK - GM_skewflx*GM_background_K
486    #endif
487  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
488       &     +0.5*(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
489  #endif  #endif
490       &     )*SlopeX(i,j)*taperFct(i,j)       &     )*SlopeX(i,j)*taperFct(i,j)
491           ENDDO           ENDDO
492          ENDDO          ENDDO
493        ENDIF         ENDIF
494  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
495    
496    #ifdef ALLOW_DIAGNOSTICS
497           IF (doDiagRediFlx) THEN
498            km1 = MAX(k-1,1)
499            DO j=1,sNy
500             DO i=1,sNx+1
501    C         store in tmp1k Kuz_Redi
502              tmp1k(i,j) = ( GM_isopycK
503    #ifdef GM_VISBECK_VARIABLE_K
504         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i-1,j,bi,bj))*0.5 _d 0
505    #endif
506         &                 )*SlopeX(i,j)*taperFct(i,j)
507             ENDDO
508            ENDDO
509            DO j=1,sNy
510             DO i=1,sNx+1
511    C-        Vertical gradients interpolated to U points
512              dTdz = (
513         &     +recip_drC(k)*
514         &       ( maskC(i-1,j,k,bi,bj)*
515         &           (theta(i-1,j,km1,bi,bj)-theta(i-1,j,k,bi,bj))
516         &        +maskC( i ,j,k,bi,bj)*
517         &           (theta( i ,j,km1,bi,bj)-theta( i ,j,k,bi,bj))
518         &       )
519         &     +recip_drC(kp1)*
520         &       ( maskC(i-1,j,kp1,bi,bj)*
521         &           (theta(i-1,j,k,bi,bj)-theta(i-1,j,kp1,bi,bj))
522         &        +maskC( i ,j,kp1,bi,bj)*
523         &           (theta( i ,j,k,bi,bj)-theta( i ,j,kp1,bi,bj))
524         &       )      ) * 0.25 _d 0
525               tmp1k(i,j) = dyG(i,j,bi,bj)*drF(k)
526         &                * _hFacW(i,j,k,bi,bj)
527         &                * tmp1k(i,j) * dTdz
528             ENDDO
529            ENDDO
530            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KuzTz', k,1,2,bi,bj,myThid)
531           ENDIF
532    #endif /* ALLOW_DIAGNOSTICS */
533    
534  C     Gradient of Sigma at V points  C     Gradient of Sigma at V points
535        DO j=1-Oly+1,sNy+Oly-1         DO j=1-Oly+1,sNy+Oly-1
536         DO i=1-Olx+1,sNx+Olx-1          DO i=1-Olx+1,sNx+Olx-1
537          SlopeX(i,j)=0.25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)           dSigmaDx(i,j)=op25*( sigmaX(i, j ,k) +sigmaX(i+1, j ,k)
538       &                    +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k) )       &                       +sigmaX(i,j-1,k) +sigmaX(i+1,j-1,k)
539       &          *_maskS(i,j,k,bi,bj)       &                      )*_maskS(i,j,k,bi,bj)
540          SlopeY(i,j)=sigmaY(i,j,k)           dSigmaDy(i,j)=sigmaY(i,j,k)
541       &          *_maskS(i,j,k,bi,bj)       &                       *_maskS(i,j,k,bi,bj)
542          dSigmaDrReal(i,j)=0.25*( sigmaR(i,j-1, k ) +sigmaR(i,j, k )           dSigmaDr(i,j)=op25*( sigmaR(i,j-1, k )+sigmaR(i,j, k )
543       &                  +maskp1*(sigmaR(i,j-1,kp1) +sigmaR(i,j,kp1)) )       &                      +(sigmaR(i,j-1,kp1)+sigmaR(i,j,kp1))*maskp1
544       &          *_maskS(i,j,k,bi,bj)       &                      )*_maskS(i,j,k,bi,bj)
545            ENDDO
546         ENDDO         ENDDO
547        ENDDO  
548    #ifdef ALLOW_AUTODIFF_TAMC
549    CADJ STORE dSigmaDx(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
550    CADJ STORE dSigmaDy(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
551    CADJ STORE dSigmaDr(:,:)   = comlev1_bibj_k, key=kkey, byte=isbyte
552    #endif /* ALLOW_AUTODIFF_TAMC */
553    
554  C     Calculate slopes for use in tensor, taper and/or clip  C     Calculate slopes for use in tensor, taper and/or clip
555        CALL GMREDI_SLOPE_LIMIT(         CALL GMREDI_SLOPE_LIMIT(
556       U             dSigmadRReal,       O             SlopeX, SlopeY,
      I             rF(K),  
      U             SlopeX, SlopeY,  
557       O             SlopeSqr, taperFct,       O             SlopeSqr, taperFct,
558         U             dSigmaDr,
559         I             dSigmaDx, dSigmaDy,
560         I             ldd97_LrhoS,rC(k),k,
561       I             bi, bj, myThid )       I             bi, bj, myThid )
562    
563    cph(
564    #ifdef ALLOW_AUTODIFF_TAMC
565    cph(
566    CADJ STORE taperfct(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
567    cph)
568    #endif /* ALLOW_AUTODIFF_TAMC */
569    cph)
570    
571  #ifdef GM_NON_UNITY_DIAGONAL  #ifdef GM_NON_UNITY_DIAGONAL
572    c      IF ( GM_nonUnitDiag ) THEN
573          DO j=1-Oly+1,sNy+Oly-1          DO j=1-Oly+1,sNy+Oly-1
574           DO i=1-Olx+1,sNx+Olx-1           DO i=1-Olx+1,sNx+Olx-1
575            Kvy(i,j,k,bi,bj) =            Kvy(i,j,k,bi,bj) =
576       &     ( GM_isopycK       &     ( GM_isopycK
577  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
578       &     +0.5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))       &     +op5*(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))
579    #endif
580         &     )
581         &     *taperFct(i,j)
582             ENDDO
583            ENDDO
584    #ifdef ALLOW_AUTODIFF_TAMC
585    # ifdef GM_EXCLUDE_CLIPPING
586    CADJ STORE Kvy(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
587    # endif
588  #endif  #endif
589       &     )*taperFct(i,j)          DO j=1-Oly+1,sNy+Oly-1
590             DO i=1-Olx+1,sNx+Olx-1
591            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 )
592           ENDDO           ENDDO
593          ENDDO          ENDDO
594    c      ENDIF
595  #endif /* GM_NON_UNITY_DIAGONAL */  #endif /* GM_NON_UNITY_DIAGONAL */
596    
597  #ifdef GM_EXTRA_DIAGONAL  #ifdef GM_EXTRA_DIAGONAL
598        IF (GM_ExtraDiag) THEN  
599    #ifdef ALLOW_AUTODIFF_TAMC
600    CADJ STORE SlopeY(:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
601    CADJ STORE taperFct(:,:)     = comlev1_bibj_k, key=kkey, byte=isbyte
602    #endif /* ALLOW_AUTODIFF_TAMC */
603           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            Kvz(i,j,k,bi,bj) =            Kvz(i,j,k,bi,bj) =
607    #ifdef ALLOW_KAPGM_CONTROL
608         &     ( GM_isopycK - GM_skewflx*kapgm(i,j,k,bi,bj)
609    #else
610       &     ( GM_isopycK - GM_skewflx*GM_background_K       &     ( GM_isopycK - GM_skewflx*GM_background_K
611    #endif
612  #ifdef GM_VISBECK_VARIABLE_K  #ifdef GM_VISBECK_VARIABLE_K
613       &     +0.5*(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
614  #endif  #endif
615       &     )*SlopeY(i,j)*taperFct(i,j)       &     )*SlopeY(i,j)*taperFct(i,j)
616           ENDDO           ENDDO
617          ENDDO          ENDDO
618        ENDIF         ENDIF
619  #endif /* GM_EXTRA_DIAGONAL */  #endif /* GM_EXTRA_DIAGONAL */
620    
621  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */  #ifdef ALLOW_DIAGNOSTICS
622           IF (doDiagRediFlx) THEN
623    c       km1 = MAX(k-1,1)
624            DO j=1,sNy+1
625             DO i=1,sNx
626    C         store in tmp1k Kvz_Redi
627              tmp1k(i,j) = ( GM_isopycK
628    #ifdef GM_VISBECK_VARIABLE_K
629         &     +(VisbeckK(i,j,bi,bj)+VisbeckK(i,j-1,bi,bj))*0.5 _d 0
630    #endif
631         &                 )*SlopeY(i,j)*taperFct(i,j)
632             ENDDO
633            ENDDO
634            DO j=1,sNy+1
635             DO i=1,sNx
636    C-        Vertical gradients interpolated to U points
637              dTdz = (
638         &     +recip_drC(k)*
639         &       ( maskC(i,j-1,k,bi,bj)*
640         &           (theta(i,j-1,km1,bi,bj)-theta(i,j-1,k,bi,bj))
641         &        +maskC(i, j ,k,bi,bj)*
642         &           (theta(i, j ,km1,bi,bj)-theta(i, j ,k,bi,bj))
643         &       )
644         &     +recip_drC(kp1)*
645         &       ( maskC(i,j-1,kp1,bi,bj)*
646         &           (theta(i,j-1,k,bi,bj)-theta(i,j-1,kp1,bi,bj))
647         &        +maskC(i, j ,kp1,bi,bj)*
648         &           (theta(i, j ,k,bi,bj)-theta(i, j ,kp1,bi,bj))
649         &       )      ) * 0.25 _d 0
650               tmp1k(i,j) = dxG(i,j,bi,bj)*drF(k)
651         &                * _hFacS(i,j,k,bi,bj)
652         &                * tmp1k(i,j) * dTdz
653             ENDDO
654            ENDDO
655            CALL DIAGNOSTICS_FILL(tmp1k, 'GM_KvzTz', k,1,2,bi,bj,myThid)
656           ENDIF
657    #endif /* ALLOW_DIAGNOSTICS */
658    
659  #ifdef ALLOW_TIMEAVE  #endif /* GM_NON_UNITY_DIAGONAL || GM_EXTRA_DIAGONAL */
 C--   Time-average  
       DO j=1-Oly+1,sNy+Oly-1  
        DO i=1-Olx+1,sNx+Olx-1  
         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  
        ENDDO  
       ENDDO  
       GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock  
 #endif /* ALLOW_TIMEAVE */  
660    
661  C-- end 2nd loop on vertical level index k  C-- end 2nd loop on vertical level index k
662        ENDDO        ENDDO
# Line 302  C-- end 2nd loop on vertical level index Line 664  C-- end 2nd loop on vertical level index
664    
665  #ifdef GM_BOLUS_ADVEC  #ifdef GM_BOLUS_ADVEC
666        IF (GM_AdvForm) THEN        IF (GM_AdvForm) THEN
667          CALL GMREDI_CALC_PSI_B(         CALL GMREDI_CALC_PSI_B(
668       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
669       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
670       I             myThid )       I             ldd97_LrhoW, ldd97_LrhoS,
671         I             myThid )
672        ENDIF        ENDIF
673  #endif  #endif
674    
675    #ifdef ALLOW_TIMEAVE
676    C--   Time-average
677          IF ( taveFreq.GT.0. ) THEN
678    
679             CALL TIMEAVE_CUMULATE( GM_Kwx_T, Kwx, Nr,
680         &                          deltaTclock, bi, bj, myThid )
681             CALL TIMEAVE_CUMULATE( GM_Kwy_T, Kwy, Nr,
682         &                          deltaTclock, bi, bj, myThid )
683             CALL TIMEAVE_CUMULATE( GM_Kwz_T, Kwz, Nr,
684         &                          deltaTclock, bi, bj, myThid )
685    #ifdef GM_VISBECK_VARIABLE_K
686           IF ( GM_Visbeck_alpha.NE.0. ) THEN
687             CALL TIMEAVE_CUMULATE( Visbeck_K_T, VisbeckK, 1,
688         &                          deltaTclock, bi, bj, myThid )
689           ENDIF
690    #endif
691    #ifdef GM_BOLUS_ADVEC
692           IF ( GM_AdvForm ) THEN
693             CALL TIMEAVE_CUMULATE( GM_PsiXtave, GM_PsiX, Nr,
694         &                          deltaTclock, bi, bj, myThid )
695             CALL TIMEAVE_CUMULATE( GM_PsiYtave, GM_PsiY, Nr,
696         &                          deltaTclock, bi, bj, myThid )
697           ENDIF
698    #endif
699           DO k=1,Nr
700             GM_TimeAve(k,bi,bj)=GM_TimeAve(k,bi,bj)+deltaTclock
701           ENDDO
702    
703          ENDIF
704    #endif /* ALLOW_TIMEAVE */
705    
706    #ifdef ALLOW_DIAGNOSTICS
707          IF ( useDiagnostics ) THEN
708            CALL GMREDI_DIAGNOSTICS_FILL(bi,bj,myThid)
709          ENDIF
710    #endif /* ALLOW_DIAGNOSTICS */
711    
712  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
713    
714        RETURN        RETURN
# Line 328  C     \================================= Line 728  C     \=================================
728    
729  C     == Global variables ==  C     == Global variables ==
730  #include "SIZE.h"  #include "SIZE.h"
 #include "GRID.h"  
 #include "DYNVARS.h"  
731  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "PARAMS.h"  
732  #include "GMREDI.h"  #include "GMREDI.h"
733    
734  C     == Routine arguments ==  C     == Routine arguments ==
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22