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

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

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

revision 1.5 by m_bates, Thu Aug 8 22:39:36 2013 UTC revision 1.22 by dfer, Thu Mar 12 15:56:02 2015 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "GMREDI_OPTIONS.h"  #include "GMREDI_OPTIONS.h"
5    
6  C     !ROUTINE: GMREDI_K3D  C     !ROUTINE: GMREDI_K3D
# Line 20  C     \ev Line 21  C     \ev
21    
22  C     == Global variables ==  C     == Global variables ==
23  #include "SIZE.h"  #include "SIZE.h"
 #include "GRID.h"  
 #include "DYNVARS.h"  
24  #include "EEPARAMS.h"  #include "EEPARAMS.h"
25  #include "PARAMS.h"  #include "PARAMS.h"
26    #include "GRID.h"
27    #include "DYNVARS.h"
28    #include "FFIELDS.h"
29  #include "GMREDI.h"  #include "GMREDI.h"
30    
31  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
# Line 47  C     === Functions ==== Line 49  C     === Functions ====
49    
50  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
51  C     == Local variables ==  C     == Local variables ==
52        INTEGER i,j,k,kk,m        INTEGER i,j,k,kk,m,kp1
53    
54  C     update_modes :: Whether to update the eigenmodes  C     update_modes :: Whether to update the eigenmodes
55        LOGICAL update_modes        LOGICAL update_modes
# Line 64  C     kLow_V :: Local version of the ind Line 66  C     kLow_V :: Local version of the ind
66  C     N2loc  :: local N**2  C     N2loc  :: local N**2
67  C     slope  :: local slope  C     slope  :: local slope
68  C     Req    :: local equatorial deformation radius (m)  C     Req    :: local equatorial deformation radius (m)
 C     Rurms  :: a local mixing length used in calculation of urms (m)  
69  C     deltaH :: local thickness of Eady integration (m)  C     deltaH :: local thickness of Eady integration (m)
70  C     g_reciprho_sq :: (gravity*recip_rhoConst)**2  C     g_reciprho_sq :: (gravity*recip_rhoConst)**2
71  C     M4loc  :: local M**4  C     M4loc  :: local M**4
# Line 74  C     sigy   :: local d(rho)/dy Line 75  C     sigy   :: local d(rho)/dy
75  C     sigz   :: local d(rho)/dz  C     sigz   :: local d(rho)/dz
76  C     hsurf  :: local surface layer depth  C     hsurf  :: local surface layer depth
77  C     small  :: a small number (to avoid floating point exceptions)  C     small  :: a small number (to avoid floating point exceptions)
78        _RL N2loc        _RL N2loc(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
79        _RL slope        _RL slope
80          _RL slopeC(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
81        _RL Req        _RL Req
82        _RL Rurms        _RL deltaH(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
       _RL deltaH  
83        _RL g_reciprho_sq        _RL g_reciprho_sq
84        _RL M4loc        _RL M4loc(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
85          _RL M4onN2(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
86        _RL maxDRhoDz        _RL maxDRhoDz
87        _RL sigx, sigy, sigz        _RL sigx, sigy, sigz
88        _RL hsurf        _RL hsurf, mskp1
89        _RL small        _RL small
90    
91  C     dfdy    :: gradient of the Coriolis paramter, df/dy, 1/(m*s)  C     dfdy    :: gradient of the Coriolis paramater, df/dy, 1/(m*s)
92  C     dfdx    :: gradient of the Coriolis paramter, df/dx, 1/(m*s)  C     dfdx    :: gradient of the Coriolis paramater, df/dx, 1/(m*s)
93  C     gradf   :: total gradient of the Coriolis paramter, SQRT(df/dx**2+df/dy**2), 1/(m*s)  C     gradf   :: gradient of the Coriolis paramater at a cell centre, 1/(m*s)
94    C     Rurms  :: a local mixing length used in calculation of urms (m)
95  C     RRhines :: The Rhines scale (m)  C     RRhines :: The Rhines scale (m)
96  C     Rmix    :: Mixing length  C     Rmix    :: Mixing length
97  C     N2      :: Square of the buoyancy frequency (1/s**2)  C     N2      :: Square of the buoyancy frequency (1/s**2)
# Line 103  C     Ubaro   :: Barotropic velocity (m/ Line 106  C     Ubaro   :: Barotropic velocity (m/
106        _RL dfdx(   1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dfdx(   1-Olx:sNx+Olx,1-Oly:sNy+Oly)
107        _RL gradf(  1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL gradf(  1-Olx:sNx+Olx,1-Oly:sNy+Oly)
108        _RL dummy(  1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dummy(  1-Olx:sNx+Olx,1-Oly:sNy+Oly)
109          _RL Rurms(  1-Olx:sNx+Olx,1-Oly:sNy+Oly)
110        _RL RRhines(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL RRhines(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
111        _RL Rmix(   1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Rmix(   1-Olx:sNx+Olx,1-Oly:sNy+Oly)
112        _RL N2(     1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL N2(     1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
# Line 111  C     Ubaro   :: Barotropic velocity (m/ Line 115  C     Ubaro   :: Barotropic velocity (m/
115        _RL N(      1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL N(      1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
116        _RL BVint(  1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL BVint(  1-Olx:sNx+Olx,1-Oly:sNy+Oly)
117        _RL Ubaro(  1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Ubaro(  1-Olx:sNx+Olx,1-Oly:sNy+Oly)
118        _RL ubar(   1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL ubar(   1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
119        _RL vbar(   1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)  
120          _RL tmpU(  1-Olx:sNx+Olx,1-Oly:sNy+Oly )
121          _RL tmpV(  1-Olx:sNx+Olx,1-Oly:sNy+Oly )
122          _RL uFldX( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr )
123          _RL vFldY( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr )
124    
125  C     Rmid     :: Rossby radius (m)  C     Rmid     :: Rossby radius (m)
126  C     KPV      :: Diffusivity (m**2/s)  C     KPV      :: Diffusivity (m**2/s)
127    C     Kdqdx    :: diffusivity multiplied by zonal PV gradient
128    C     Kdqdy    :: diffusivity multiplied by meridional PV gradient
129  C     SlopeX   :: isopycnal slope in x direction  C     SlopeX   :: isopycnal slope in x direction
130  C     SlopeY   :: isopycnal slope in y direction  C     SlopeY   :: isopycnal slope in y direction
131  C     dSigmaDx :: sigmaX averaged onto tracer grid  C     dSigmaDx :: sigmaX averaged onto tracer grid
# Line 130  C     coriV    :: As for cori, but, at V Line 140  C     coriV    :: As for cori, but, at V
140  C     surfkz   :: Depth of surface layer (in r units)  C     surfkz   :: Depth of surface layer (in r units)
141        _RL Rmid(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Rmid(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
142        _RL KPV(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KPV(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
143          _RL Kdqdy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
144          _RL Kdqdx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
145        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
146        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
147        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
# Line 143  C     surfkz   :: Depth of surface layer Line 155  C     surfkz   :: Depth of surface layer
155        _RL fCoriV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL fCoriV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
156        _RL surfkz(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL surfkz(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
157    
158    C     centreX,centreY :: used for calculating averages at centre of cell
159    C     numerator,denominator :: of the renormalisation factor
160    C     uInt      :: column integral of u velocity (sum u*dz)
161    C     vInt      :: column integral of v velocity (sum v*dz)
162    C     KdqdxInt  :: column integral of K*dqdx (sum K*dqdx*dz)
163    C     KdqdyInt  :: column integral of K*dqdy (sum K*dqdy*dz)
164    C     uKdqdyInt :: column integral of u*K*dqdy (sum u*K*dqdy*dz)
165    C     vKdqdxInt :: column integral of v*K*dqdx (sum v*K*dqdx*dz)
166    C     uXiyInt   :: column integral of u*Xiy (sum u*Xiy*dz)
167    C     vXixInt   :: column integral of v*Xix (sum v*Xix*dz)
168    C     Renorm    :: renormalisation factor at the centre of a cell
169    C     RenormU   :: renormalisation factor at the western face of a cell
170    C     RenormV   :: renormalisation factor at the southern face of a cell
171          _RL centreX, centreY
172          _RL numerator, denominator
173          _RL uInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
174          _RL vInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
175          _RL KdqdxInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
176          _RL KdqdyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
177          _RL uKdqdyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
178          _RL vKdqdxInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
179          _RL uXiyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
180          _RL vXixInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
181          _RL Renorm(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
182          _RL RenormU(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
183          _RL RenormV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
184    
185  C     gradqx   :: Potential vorticity gradient in x direction  C     gradqx   :: Potential vorticity gradient in x direction
186  C     gradqy   :: Potential vorticity gradient in y direction  C     gradqy   :: Potential vorticity gradient in y direction
187  C     XimX     :: Vertical integral of phi_m*K*gradqx  C     XimX     :: Vertical integral of phi_m*K*gradqx
# Line 174  C     psistar :: eddy induced streamfunc Line 213  C     psistar :: eddy induced streamfunc
213        _RL psistar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)        _RL psistar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
214  #endif  #endif
215    
   
216  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
217    
218  C     ======================================  C     ======================================
# Line 201  C     ================================== Line 239  C     ==================================
239         ENDDO         ENDDO
240        ENDDO        ENDDO
241    
242  C     Dummy values for the edges  C     Dummy values for the edges. This does not affect the results
243  C     This avoids weirdness in gmredi_calc_eigs  C     but avoids problems when solving for the eigenvalues.
244        i=1-Olx        i=1-Olx
245        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
246         kLow_U(i,j) = kLow_C(i,j)         kLow_U(i,j) = 0
247        ENDDO        ENDDO
248        j=1-Oly        j=1-Oly
249        DO i=1-Olx,sNx+Olx        DO i=1-Olx,sNx+Olx
250         kLow_V(i,j) = kLow_C(i,j)         kLow_V(i,j) = 0
251        ENDDO        ENDDO
252    
253        g_reciprho_sq = (gravity*recip_rhoConst)**2        g_reciprho_sq = (gravity*recip_rhoConst)**2
# Line 223  C     Gradient of Coriolis Line 261  C     Gradient of Coriolis
261         ENDDO         ENDDO
262        ENDDO        ENDDO
263    
264  C     Coriolis at C points enforcing a minimum value so  C     Coriolis at C points enforcing a minimum value so
265  C     that it is defined at the equator  C     that it is defined at the equator
266        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
267         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
# Line 232  C     that it is defined at the equator Line 270  C     that it is defined at the equator
270         ENDDO         ENDDO
271        ENDDO        ENDDO
272  C     Coriolis at U and V points  C     Coriolis at U and V points
273        DO j=1-Oly+1,sNy+Oly        DO j=1-Oly,sNy+Oly
274         DO i=1-Olx+1,sNx+Olx         DO i=1-Olx+1,sNx+Olx
275  C       Limited so that the inverse is defined at the equator  C       Limited so that the inverse is defined at the equator
276          coriU(i,j) = op5*( cori(i,j)+cori(i-1,j) )          coriU(i,j) = op5*( cori(i,j)+cori(i-1,j) )
277          coriU(i,j) = SIGN( MAX( ABS(coriU(i,j)),GM_K3D_minCori ),          coriU(i,j) = SIGN( MAX( ABS(coriU(i,j)),GM_K3D_minCori ),
278       &       coriU(i,j) )       &       coriU(i,j) )
279    
280    C       Not limited
281            fCoriU(i,j) = op5*( fCori(i,j,bi,bj)+fCori(i-1,j,bi,bj) )
282           ENDDO
283          ENDDO
284          DO j=1-Oly+1,sNy+Oly
285           DO i=1-Olx,sNx+Olx
286    C       Limited so that the inverse is defined at the equator
287          coriV(i,j) = op5*( cori(i,j)+cori(i,j-1) )          coriV(i,j) = op5*( cori(i,j)+cori(i,j-1) )
288          coriV(i,j) = SIGN( MAX( ABS(coriV(i,j)),GM_K3D_minCori ),          coriV(i,j) = SIGN( MAX( ABS(coriV(i,j)),GM_K3D_minCori ),
289       &       coriV(i,j) )       &       coriV(i,j) )
290    
291  C       Not limited  C       Not limited
         fCoriU(i,j) = op5*( fCori(i,j,bi,bj)+fCori(i-1,j,bi,bj) )  
292          fCoriV(i,j) = op5*( fCori(i,j,bi,bj)+fCori(i,j-1,bi,bj) )          fCoriV(i,j) = op5*( fCori(i,j,bi,bj)+fCori(i,j-1,bi,bj) )
293         ENDDO         ENDDO
294        ENDDO        ENDDO
295    C     Computing beta
296          IF ( selectCoriMap.EQ.1 ) THEN
297           DO j=1-Oly,sNy+Oly
298            DO i=1-Olx,sNx+Olx
299             gradf(i,j) =  beta
300            ENDDO
301           ENDDO
302          ELSEIF ( selectCoriMap.EQ.2 ) THEN
303           DO j=1-Oly,sNy+Oly
304            DO i=1-Olx,sNx+Olx
305             gradf(i,j) = recip_rSphere*fCoriCos(i,j,bi,bj)
306            ENDDO
307           ENDDO
308          ENDIF
309    C     Some dummy values at the edges
310          i=1-Olx
311          DO j=1-Oly,sNy+Oly
312           coriU(i,j)=cori(i,j)
313           fCoriU(i,j)=fCori(i,j,bi,bj)
314          ENDDO
315          j=1-Oly
316          DO i=1-Olx,sNx+Olx
317           coriV(i,j)=cori(i,j)
318           fCoriV(i,j)=fCori(i,j,bi,bj)
319          ENDDO
320    
321    C     Zeroing some cumulative fields
322        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
323         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
324          gradf(i,j) = SQRT( dfdx(i,j)**2 + dfdy(i,j)**2 )          eady(i,j)   = zeroRL
325            BVint(i,j)  = zeroRL
326            Ubaro(i,j)  = zeroRL
327            deltaH(i,j) = zeroRL
328           ENDDO
329          ENDDO
330          DO k=1,Nr
331           DO j=1-Oly,sNy+Oly
332            DO i=1-Olx,sNx+Olx
333             slopeC(i,j,k)=zeroRL
334            ENDDO
335         ENDDO         ENDDO
336        ENDDO        ENDDO
337    
338  C     Zeroing some cumulative fields  C     initialise remaining 2d variables
339        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
340         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
341          eady(i,j)  = zeroRL          dfdy(i,j)=zeroRL
342          BVint(i,j) = zeroRL          dfdy(i,j)=zeroRL
343          Ubaro(i,j) = zeroRL          Rurms(i,j)=zeroRL
344            RRhines(i,j)=zeroRL
345            Rmix(i,j)=zeroRL
346           ENDDO
347          ENDDO
348    C     initialise remaining 3d variables
349          DO k=1,Nr
350           DO j=1-Oly,sNy+Oly
351            DO i=1-Olx,sNx+Olx
352             N2loc(i,j,k)=GM_K3D_minN2
353             N2W(i,j,k) = GM_K3D_minN2
354             N2S(i,j,k) = GM_K3D_minN2
355             M4loc(i,j,k)=zeroRL
356             M4onN2(i,j,k)=zeroRL
357             urms(i,j,k)=zeroRL
358             SlopeX(i,j,k)=zeroRL
359             SlopeY(i,j,k)=zeroRL
360             dSigmaDx(i,j,k)=zeroRL
361             dSigmaDy(i,j,k)=zeroRL
362             gradqx(i,j,k)=zeroRL
363             gradqy(i,j,k)=zeroRL
364            ENDDO
365         ENDDO         ENDDO
366        ENDDO        ENDDO
367    
368  C     Find the zonal velocity at the cell centre  C     Find the zonal velocity at the cell centre
369  C     The logicals here are, in order: 1/ go from grid to north/east directions  #ifdef ALLOW_EDDYPSI
370  C     2/ go from C to A grid and 3/ apply the mask        IF (GM_InMomAsStress) THEN
371        CALL rotate_uv2en_rl(uVel, vVel, ubar, vbar, .TRUE., .TRUE.,          DO k=1,Nr
372       &     .TRUE.,Nr,mythid)           DO i = 1-olx,snx+olx
373              DO j = 1-oly,sny+oly
374               uFldX(i,j,k) = uEulerMean(i,j,k,bi,bj)
375               vFldY(i,j,k) = vEulerMean(i,j,k,bi,bj)
376              ENDDO
377             ENDDO
378            ENDDO
379          ELSE
380    #endif
381            DO k=1,Nr
382             DO i = 1-olx,snx+olx
383              DO j = 1-oly,sny+oly
384               uFldX(i,j,k) = uVel(i,j,k,bi,bj)
385               vFldY(i,j,k) = vVel(i,j,k,bi,bj)
386              ENDDO
387             ENDDO
388            ENDDO
389    #ifdef ALLOW_EDDYPSI
390          ENDIF
391    #endif
392    
393    C     The following comes from rotate_uv2en_rl
394    C     This code does two things:
395    C     1) go from C grid velocity points to A grid velocity points
396    C     2) go from model grid directions to east/west directions
397          DO k = 1,Nr
398    
399           DO i = 1-Olx,sNx+Olx
400            j=sNy+Oly
401            tmpU(i,j)=zeroRL
402            tmpV(i,j)=zeroRL
403           ENDDO
404           DO j = 1-Oly,sNy+Oly-1
405            i=sNx+Olx
406            tmpU(i,j)=zeroRL
407            tmpV(i,j)=zeroRL
408            DO i = 1-Olx,sNx+Olx-1
409             tmpU(i,j) = 0.5 _d 0
410         &        *( uFldX(i+1,j,k) + uFldX(i,j,k) )
411             tmpV(i,j) = 0.5 _d 0
412         &        *( vFldY(i,j+1,k) + vFldY(i,j,k) )
413    
414             tmpU(i,j) = tmpU(i,j) * maskC(i,j,k,bi,bj)
415             tmpV(i,j) = tmpV(i,j) * maskC(i,j,k,bi,bj)
416            ENDDO
417           ENDDO
418    
419           DO j = 1-oly,sny+oly
420            DO i = 1-olx,snx+olx
421            ENDDO
422           ENDDO
423    
424    C     rotation
425           DO j = 1-oly,sny+oly
426            DO i = 1-olx,snx+olx
427             ubar(i,j,k) =
428         &        angleCosC(i,j,bi,bj)*tmpU(i,j)
429         &        -angleSinC(i,j,bi,bj)*tmpV(i,j)
430            ENDDO
431           ENDDO
432          ENDDO
433    
434  C     Square of the buoyancy frequency at the top of a grid cell  C     Square of the buoyancy frequency at the top of a grid cell
435    C     Enforce a minimum N2
436    C     Mask N2, so it is zero at bottom
437        DO k=2,Nr        DO k=2,Nr
438         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
439          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
440           N2(i,j,k) = -gravity*recip_rhoConst*sigmaR(i,j,k)           N2(i,j,k) = -gravity*recip_rhoConst*sigmaR(i,j,k)
441             N2(i,j,k) = MAX(N2(i,j,k),GM_K3D_minN2)*maskC(i,j,k,bi,bj)
442             N(i,j,k)  = SQRT(N2(i,j,k))
443          ENDDO          ENDDO
444         ENDDO         ENDDO
445        ENDDO        ENDDO
446  C     N2(k=1) is always zero  C     N2(k=1) is always zero
       k=1  
447        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
448         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
449          N2(i,j,k) = 0.0          N2(i,j,1) = zeroRL
450          N(i,j,k)  = 0.0          N(i,j,1)  = zeroRL
        ENDDO  
       ENDDO  
 C     Enforce a minimum N2  
       DO k=2,Nr  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          IF (N2(i,j,k).LT.GM_K3D_minN2) N2(i,j,k)=GM_K3D_minN2  
          N(i,j,k) = SQRT(N2(i,j,k))  
         ENDDO  
451         ENDDO         ENDDO
452        ENDDO        ENDDO
453  C     Calculate the minimum drho/dz  C     Calculate the minimum drho/dz
454        maxDRhoDz = -rhoConst*GM_K3D_minN2/gravity        maxDRhoDz = -rhoConst*GM_K3D_minN2/gravity
455    
456  C     Calculate the barotropic velocity by vertically integrating  C     Calculate the barotropic velocity by vertically integrating
457  C     and the dividing by the depth of the water column  C     and the dividing by the depth of the water column
458  C     Note that Ubaro is on the U grid.  C     Note that Ubaro is at the C-point.
459        DO k=1,Nr        DO k=1,Nr
460         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
461          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
462           Ubaro(i,j) = Ubaro(i,j) +           Ubaro(i,j) = Ubaro(i,j) +
463       &        maskW(i,j,k,bi,bj)*drF(k)*hfacC(i,j,k,bi,bj)       &        drF(k)*hfacC(i,j,k,bi,bj)*ubar(i,j,k)
      &                *ubar(i,j,k,bi,bj)  
464          ENDDO          ENDDO
465         ENDDO         ENDDO
466        ENDDO        ENDDO
# Line 319  C         The minus sign is because r_Lo Line 474  C         The minus sign is because r_Lo
474        ENDDO        ENDDO
475    
476  C     Integrate the buoyancy frequency vertically using the trapezoidal method.  C     Integrate the buoyancy frequency vertically using the trapezoidal method.
477    C     Assume that N(z=-H)=0
478        DO k=1,Nr        DO k=1,Nr
479           kp1 = min(k+1,Nr)
480           mskp1 = oneRL
481           IF ( k.EQ.Nr ) mskp1 = zeroRL
482         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
483          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
          IF (k.LT.kLow_C(i,j)) THEN  
484             BVint(i,j) = BVint(i,j) + hFacC(i,j,k,bi,bj)*drF(k)             BVint(i,j) = BVint(i,j) + hFacC(i,j,k,bi,bj)*drF(k)
485       &                               *(N(i,j,k)+N(i,j,k+1))       &                         *op5*(N(i,j,k)+mskp1*N(i,j,kp1))
          ELSEIF (k.EQ.kLow_C(i,j)) THEN  
 C          Assume that N(z=-H)=0  
            BVint(i,j) = BVint(i,j) + hFacC(i,j,k,bi,bj)*drF(k)*N(i,j,k)  
          ENDIF  
486          ENDDO          ENDDO
487         ENDDO         ENDDO
488        ENDDO        ENDDO
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         BVint(i,j) = op5*BVint(i,j)  
        ENDDO  
       ENDDO  
489    
490  C     Calculate the eigenvalues and eigenvectors  C     Calculate the eigenvalues and eigenvectors
491        IF (update_modes) THEN        IF (update_modes) THEN
# Line 357  C       Calculate the Rossby Radius Line 506  C       Calculate the Rossby Radius
506        ENDIF        ENDIF
507    
508  C     Average dsigma/dx and dsigma/dy onto the centre points  C     Average dsigma/dx and dsigma/dy onto the centre points
509          
510        DO k=1,Nr        DO k=1,Nr
511         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
512          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
# Line 372  C     Calculate the Eady growth rate Line 521  C     Calculate the Eady growth rate
521  C     ===============================  C     ===============================
522        DO k=1,Nr        DO k=1,Nr
523    
524           kp1 = min(k+1,Nr)
525           mskp1 = oneRL
526           IF ( k.EQ.Nr ) mskp1 = zeroRL
527    
528           DO j=1-Oly,sNy+Oly-1
529            DO i=1-Olx,sNx+Olx-1
530             M4loc(i,j,k) = g_reciprho_sq*( dSigmaDx(i,j,k)**2
531         &                                 +dSigmaDy(i,j,k)**2 )
532             N2loc(i,j,k) = op5*(N2(i,j,k)+mskp1*N2(i,j,kp1))
533            ENDDO
534           ENDDO
535  C      The bottom of the grid cell is shallower than the top  C      The bottom of the grid cell is shallower than the top
536  C      integration level, so, advance the depth.  C      integration level, so, advance the depth.
537         IF (-rF(k+1).LE. GM_K3D_EadyMinDepth) CYCLE         IF (-rF(k+1) .LE. GM_K3D_EadyMinDepth) CYCLE
538    
539  C      Don't bother going any deeper since the top of the  C      Do not bother going any deeper since the top of the
540  C      cell is deeper than the bottom integration level  C      cell is deeper than the bottom integration level
541         IF (-rF(k).GE.GM_K3D_EadyMaxDepth) EXIT         IF (-rF(k).GE.GM_K3D_EadyMaxDepth) EXIT
542    
543  C      We are in the integration depth range  C      We are in the integration depth range
544         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
545          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
546           IF (kLow_C(i,j).GE.k) THEN           IF ( (kLow_C(i,j).GE.k) .AND.
547             IF (k.NE.kLow_C(i,j)) THEN       &        (-hMixLayer(i,j,bi,bj).LE.-rC(k)) ) THEN
              N2loc = op5*(N2(i,j,k)+N2(i,j,k+1))  
            ELSE  
              N2loc = op5*N2(i,j,k)  
            ENDIF  
            M4loc = g_reciprho_sq*( dSigmaDx(i,j,k)**2  
      &                            +dSigmaDy(i,j,k)**2 )  
            slope = SQRT(SQRT(M4loc)/N2loc)  
548    
549               slopeC(i,j,k) = SQRT(M4loc(i,j,k))/N2loc(i,j,k)
550  C          Limit the slope.  Note, this is not all the Eady calculations.  C          Limit the slope.  Note, this is not all the Eady calculations.
551             IF (slope.LE.GM_K3D_maxSlope) THEN             IF (slopeC(i,j,k).LE.GM_maxSlope) THEN
552               eady(i,j) = eady(i,j)               M4onN2(i,j,k) = M4loc(i,j,k)/N2loc(i,j,k)
      &            + hfacC(i,j,k,bi,bj)*drF(k)*M4loc/(N2loc)  
553             ELSE             ELSE
554               eady(i,j) = eady(i,j)               slopeC(i,j,k) = GM_maxslope
555       &            + hfacC(i,j,k,bi,bj)*drF(k)*SQRT(M4loc)               M4onN2(i,j,k) = SQRT(M4loc(i,j,k))*GM_maxslope
      &              *GM_K3D_maxSlope*GM_K3D_maxSlope  
556             ENDIF             ENDIF
557               eady(i,j)   = eady(i,j)
558         &                   + hfacC(i,j,k,bi,bj)*drF(k)*M4onN2(i,j,k)
559               deltaH(i,j) = deltaH(i,j) + drF(k)
560           ENDIF           ENDIF
561          ENDDO          ENDDO
562         ENDDO         ENDDO
# Line 409  C          Limit the slope.  Note, this Line 564  C          Limit the slope.  Note, this
564    
565        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
566         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
567  C       If the minimum depth for the integration is deeper than ocean  C       If the minimum depth for the integration is deeper than the ocean
568  C       bottom then give the eady growth rate a dummy, non-zero value  C       bottom OR the mixed layer is deeper than the maximum depth of
569  C       to avoid floating point exceptions.  These points are taken care  C       integration, we set the Eady growth rate to something small
570  C       of by setting K3D=GM_K3D_smallK later.  C       to avoid floating point exceptions.
571          IF (kLow_C(i,j).NE.0  C       Later, these areas will be given a small diffusivity.
572       &       .AND. -r_Low(i,j,bi,bj).LT.GM_K3D_EadyMinDepth) THEN          IF (deltaH(i,j).EQ.zeroRL) THEN
573            eady(i,j) = small            eady(i,j) = small
574    
575  C       Otherwise, multiply eady by the various constants to get the  C       Otherwise, divide over the integration and take the square root
576  C       growth rate.  C       to actually find the Eady growth rate.
577          ELSE          ELSE
578            deltaH = MIN(-r_low(i,j,bi,bj),GM_K3D_EadyMaxDepth)            eady(i,j) = SQRT(eady(i,j)/deltaH(i,j))
579            deltaH = deltaH - GM_K3D_EadyMinDepth  
           eady(i,j) = SQRT(eady(i,j)/deltaH)  
             
580          ENDIF          ENDIF
581    
582         ENDDO         ENDDO
# Line 435  C     ================================== Line 588  C     ==================================
588        DO j=1-Oly+1,sNy+Oly        DO j=1-Oly+1,sNy+Oly
589         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
590  C       Calculate the Visbeck velocity  C       Calculate the Visbeck velocity
591          Rurms = MIN(Rdef(i,j,bi,bj),GM_K3D_maxLurms)          Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_K3D_Rmax)
592          urms(i,j,1) = GM_K3D_Lambda*eady(i,j)*Rurms          urms(i,j,1) = GM_K3D_Lambda*eady(i,j)*Rurms(i,j)
593  C       Set the bottom urms to zero  C       Set the bottom urms to zero
594          k=kLow_C(i,j)          k=kLow_C(i,j)
595          IF (k.GT.0) urms(i,j,k) = 0.0          IF (k.GT.0) urms(i,j,k) = 0.0
# Line 446  C       Calculate the Rhines scale Line 599  C       Calculate the Rhines scale
599    
600  C       Calculate the estimated length scale  C       Calculate the estimated length scale
601          Rmix(i,j) = MIN(Rdef(i,j,bi,bj), RRhines(i,j))          Rmix(i,j) = MIN(Rdef(i,j,bi,bj), RRhines(i,j))
602            Rmix(i,j) = MAX(Rmix(i,j),GM_K3D_Rmin)
603    
604  C       Calculate the Doppler shifted long Rossby wave speed  C       Calculate the Doppler shifted long Rossby wave speed
605  C       Ubaro is on the U grid so we must average onto the M grid.  C       Ubaro is at the C-point.
606          cDopp(i,j) = op5*( Ubaro(i,j)+Ubaro(i+1,j) )          cDopp(i,j) = Ubaro(i,j)
607       &                - gradf(i,j)*Rdef(i,j,bi,bj)*Rdef(i,j,bi,bj)       &                - gradf(i,j)*Rdef(i,j,bi,bj)*Rdef(i,j,bi,bj)
608  C       Limit the wave speed to the namelist variable GM_K3D_maxC  C       Limit the wave speed to the namelist variable GM_K3D_maxC
609          IF (ABS(cDopp(i,j)).GT.GM_K3D_maxC) THEN          IF (ABS(cDopp(i,j)).GT.GM_K3D_maxC) THEN
# Line 469  C     Calculate the diffusivity (on the Line 623  C     Calculate the diffusivity (on the
623         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
624          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
625           IF (k.LE.kLow_C(i,j)) THEN           IF (k.LE.kLow_C(i,j)) THEN
626             IF (-r_Low(i,j,bi,bj).LT.GM_K3D_EadyMinDepth) THEN             IF (deltaH(i,j).EQ.zeroRL) THEN
627               K3D(i,j,k,bi,bj) = GM_K3D_smallK               K3D(i,j,k,bi,bj) = GM_K3D_smallK
628             ELSE             ELSE
629               IF (urms(i,j,k).EQ.0.0) THEN               IF (urms(i,j,k).EQ.0.0) THEN
630                 K3D(i,j,k,bi,bj) = GM_K3D_smallK                 K3D(i,j,k,bi,bj) = GM_K3D_smallK
631               ELSE               ELSE
632                 umc(i,j,k)  = ubar(i,j,k,bi,bj) - cDopp(i,j)                umc(i,j,k) =ubar(i,j,k) - cDopp(i,j)
633                 supp(i,j,k) = 1/( 1 + 4*umc(i,j,k)**2/urms(i,j,1)**2 )                supp(i,j,k)=1./(1.+GM_K3D_b1*umc(i,j,k)**2/urms(i,j,1)**2)
634                 K3D(i,j,k,bi,bj) = GM_K3D_gamma*urms(i,j,k)  C              2*Rmix gives the diameter
635       &                            *Rmix(i,j)*supp(i,j,k)                K3D(i,j,k,bi,bj) = GM_K3D_gamma*urms(i,j,k)
636         &                           *2.*Rmix(i,j)*supp(i,j,k)
637               ENDIF               ENDIF
638    
639  C            Enforce lower and upper bounds on the diffusivity  C            Enforce lower and upper bounds on the diffusivity
640               IF (K3D(i,j,k,bi,bj).LT.GM_K3D_smallK)               K3D(i,j,k,bi,bj) = MIN(K3D(i,j,k,bi,bj),GM_maxK3D)
641       &            K3D(i,j,k,bi,bj) = GM_K3D_smallK               K3D(i,j,k,bi,bj) = MAX(K3D(i,j,k,bi,bj),GM_K3D_smallK)
              IF (K3D(i,j,k,bi,bj).GT.GM_maxK3D)  
      &            K3D(i,j,k,bi,bj) = GM_maxK3D  
642             ENDIF             ENDIF
643           ENDIF           ENDIF
644          ENDDO          ENDDO
# Line 496  C     ================================== Line 649  C     ==================================
649  C     Find the PV gradient  C     Find the PV gradient
650  C     ======================================  C     ======================================
651  C     Calculate the surface layer thickness.  C     Calculate the surface layer thickness.
652  C     Use hMixLayer (calculated in model/src/calc_oce_mxlayer)  C     Use hMixLayer (calculated in model/src/calc_oce_mxlayer)
653  C     for the mixed layer depth.  C     for the mixed layer depth.
654    
655  C     Enforce a minimum surface layer depth  C     Enforce a minimum surface layer depth
# Line 511  C     Enforce a minimum surface layer de Line 664  C     Enforce a minimum surface layer de
664        DO k=1,Nr        DO k=1,Nr
665         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
666          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
667           IF (rF(k).GT.surfkz(i,j) .AND. surfkz(i,j).GE.rF(k+1))           IF (rF(k).GT.surfkz(i,j) .AND. surfkz(i,j).GE.rF(k+1))
668       &        surfk(i,j) = k       &        surfk(i,j) = k
669          ENDDO          ENDDO
670         ENDDO         ENDDO
# Line 538  C          Calculate the zonal slope at Line 691  C          Calculate the zonal slope at
691             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i-1,j,k)), maxDRhoDz )             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i-1,j,k)), maxDRhoDz )
692             sigx  = op5*( sigmaX(i,j,k)+sigmaX(i,j,k-1) )             sigx  = op5*( sigmaX(i,j,k)+sigmaX(i,j,k-1) )
693             slope = sigx/sigz             slope = sigx/sigz
 C           IF(ABS(slope).GT.GM_K3D_maxSlope)  
 C     &          slope = SIGN(GM_K3D_maxSlope,slope)  
694             IF(ABS(slope).GT.GM_maxSlope)             IF(ABS(slope).GT.GM_maxSlope)
695       &          slope = SIGN(GM_maxSlope,slope)       &          slope = SIGN(GM_maxSlope,slope)
696             SlopeX(i,j,k)=-maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)*slope             SlopeX(i,j,k)=-maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)*slope
697              
698  C          Calculate the meridional slope at the southern cell face (V grid)  C          Calculate the meridional slope at the southern cell face (V grid)
699             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i,j-1,k)), maxDRhoDz )             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i,j-1,k)), maxDRhoDz )
700             sigy  = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )             sigy  = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )
701             slope = sigy/sigz             slope = sigy/sigz
 C           IF(ABS(slope).GT.GM_K3D_maxSlope)  
 C     &          slope = SIGN(GM_K3D_maxSlope,slope)  
702             IF(ABS(slope).GT.GM_maxSlope)             IF(ABS(slope).GT.GM_maxSlope)
703       &          slope = SIGN(GM_maxSlope,slope)       &          slope = SIGN(GM_maxSlope,slope)
704             SlopeY(i,j,k)=-maskS(i,j,k-1,bi,bj)*maskS(i,j,k,bi,bj)*slope             SlopeY(i,j,k)=-maskS(i,j,k-1,bi,bj)*maskS(i,j,k,bi,bj)*slope
# Line 558  C     &          slope = SIGN(GM_K3D_max Line 707  C     &          slope = SIGN(GM_K3D_max
707         ENDDO         ENDDO
708        ENDDO        ENDDO
709    
710  C     Calculate the thickness flux  C     Calculate the thickness flux and diffusivity.  These may be altered later
711    C     depending on namelist options.
712  C     Enforce a zero slope bottom boundary condition for the bottom most cells (k=Nr)  C     Enforce a zero slope bottom boundary condition for the bottom most cells (k=Nr)
713        k=Nr        k=Nr
714        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
# Line 569  C          Zonal thickness flux at the w Line 719  C          Zonal thickness flux at the w
719  C          Meridional thickness flux at the southern cell face  C          Meridional thickness flux at the southern cell face
720             tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)             tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)
721       &                     *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)       &                     *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
722    
723    C          Use the interior diffusivity. Note: if we are using a
724    C          constant diffusivity KPV is overwritten later
725               KPV(i,j,k) = K3D(i,j,k,bi,bj)
726    
727         ENDDO         ENDDO
728        ENDDO        ENDDO
729    
730  C     Calculate the thickness flux for other cells (k<Nr)  C     Calculate the thickness flux and diffusivity and for other cells (k<Nr)
 C     SlopeX and SlopeY are zero in dry cells, so, the bottom boundary  
 C     condition that the slope is zero is taken care of.  
 C     We still need to give special treatment for the surface layer however.  
731        DO k=Nr-1,1,-1        DO k=Nr-1,1,-1
732         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly
733          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx
734           IF(k.LE.surfk(i,j) .AND. .NOT. GM_K3D_likeGM) THEN  C        Zonal thickness flux at the western cell face
735  C          We are in the surface layer, so set the thickness flux           tfluxX(i,j,k)=-fCoriU(i,j)*(SlopeX(i,j,k)-SlopeX(i,j,k+1))
736  C          based on the average slope over the surface layer       &        *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
737  C          If we are on the edge of a "cliff" the surface layer at the       &        *maskW(i,j,k,bi,bj)
 C          centre of the grid point could be deeper than the U or V point.    
 C          So, we ensure that we always take a sensible slope.  
            IF(kLow_U(i,j).LT.surfk(i,j)) THEN  
              kk=kLow_U(i,j)  
              hsurf = -rLowW(i,j,bi,bj)  
            ELSE  
              kk=surfk(i,j)  
              hsurf = -surfkz(i,j)  
            ENDIF  
            IF(kk.GT.0) THEN  
              IF(kk.EQ.Nr) THEN  
                tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)  
      &              *SlopeX(i,j,kk)/hsurf  
              ELSE  
                tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)  
      &              *( SlopeX(i,j,kk)-SlopeX(i,j,kk+1) )/hsurf  
              ENDIF  
            ELSE  
              tfluxX(i,j,k) = zeroRL  
            ENDIF  
   
            IF(kLow_V(i,j).LT.surfk(i,j)) THEN  
              kk=kLow_V(i,j)  
              hsurf = -rLowS(i,j,bi,bj)  
            ELSE  
              kk=surfk(i,j)  
              hsurf = -surfkz(i,j)  
            ENDIF  
            IF(kk.GT.0) THEN  
              IF(kk.EQ.Nr) THEN  
                tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)  
      &              *SlopeY(i,j,kk)/hsurf  
              ELSE  
                tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)  
      &              *( SlopeY(i,j,kk)-SlopeY(i,j,kk+1) )/hsurf  
              ENDIF  
            ELSE  
              tfluxY(i,j,k) = zeroRL  
            ENDIF  
738    
739           ELSE  C        Meridional thickness flux at the southern cell face
740  C          We are not in the surface layer, so calculate the thickness           tfluxY(i,j,k)=-fCoriV(i,j)*(SlopeY(i,j,k)-SlopeY(i,j,k+1))
741  C          flux based on the local isopyncal slope       &        *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
742         &        *maskS(i,j,k,bi,bj)
743    
744  C          Zonal thickness flux at the western cell face  C        Use the interior diffusivity. Note: if we are using a
745             tfluxX(i,j,k)=-fCoriU(i,j)*(SlopeX(i,j,k)-SlopeX(i,j,k+1))  C        constant diffusivity KPV is overwritten later
746       &                    *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)           KPV(i,j,k) = K3D(i,j,k,bi,bj)
      &                    *maskW(i,j,k,bi,bj)  
   
 C          Meridional thickness flux at the southern cell face  
            tfluxY(i,j,k)=-fCoriV(i,j)*(SlopeY(i,j,k)-SlopeY(i,j,k+1))  
      &                    *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)  
      &                    *maskS(i,j,k,bi,bj)  
          ENDIF  
747          ENDDO          ENDDO
748         ENDDO         ENDDO
749        ENDDO        ENDDO
750    
751    C     Special treatment for the surface layer (if necessary), which overwrites
752    C     values in the previous loops.
753          IF (GM_K3D_ThickSheet .OR. GM_K3D_surfK) THEN
754            DO k=Nr-1,1,-1
755             DO j=1-Oly,sNy+Oly
756              DO i=1-Olx,sNx+Olx
757               IF(k.LE.surfk(i,j)) THEN
758    C            We are in the surface layer.  Change the thickness flux
759    C            and diffusivity as necessary.
760    
761                 IF (GM_K3D_ThickSheet) THEN
762    C              We are in the surface layer, so set the thickness flux
763    C              based on the average slope over the surface layer
764    C              If we are on the edge of a "cliff" the surface layer at the
765    C              centre of the grid point could be deeper than the U or V point.
766    C              So, we ensure that we always take a sensible slope.
767                   IF(kLow_U(i,j).LT.surfk(i,j)) THEN
768                     kk=kLow_U(i,j)
769                     hsurf = -rLowW(i,j,bi,bj)
770                   ELSE
771                     kk=surfk(i,j)
772                     hsurf = -surfkz(i,j)
773                   ENDIF
774                   IF(kk.GT.0) THEN
775                     IF(kk.EQ.Nr) THEN
776                       tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
777         &                  *SlopeX(i,j,kk)/hsurf
778                     ELSE
779                       tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
780         &                  *( SlopeX(i,j,kk)-SlopeX(i,j,kk+1) )/hsurf
781                     ENDIF
782                   ELSE
783                     tfluxX(i,j,k) = zeroRL
784                   ENDIF
785    
786                   IF(kLow_V(i,j).LT.surfk(i,j)) THEN
787                     kk=kLow_V(i,j)
788                     hsurf = -rLowS(i,j,bi,bj)
789                   ELSE
790                     kk=surfk(i,j)
791                     hsurf = -surfkz(i,j)
792                   ENDIF
793                   IF(kk.GT.0) THEN
794                     IF(kk.EQ.Nr) THEN
795                       tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
796         &                  *SlopeY(i,j,kk)/hsurf
797                     ELSE
798                       tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
799         &                  *( SlopeY(i,j,kk)-SlopeY(i,j,kk+1) )/hsurf
800                     ENDIF
801                   ELSE
802                     tfluxY(i,j,k) = zeroRL
803                   ENDIF
804                 ENDIF
805    
806                 IF (GM_K3D_surfK) THEN
807    C              Use a constant K in the surface layer.
808                   KPV(i,j,k) = GM_K3D_constK
809                 ENDIF
810               ENDIF
811              ENDDO
812             ENDDO
813            ENDDO
814          ENDIF
815    
816  C     Calculate gradq  C     Calculate gradq
817        IF (GM_K3D_likeGM .OR. GM_K3D_beta_eq_0) THEN        IF (GM_K3D_beta_eq_0) THEN
818  C       Ignore beta in the calculation of grad(q)  C       Ignore beta in the calculation of grad(q)
819          DO k=1,Nr          DO k=1,Nr
820           DO j=1-Oly+1,sNy+Oly           DO j=1-Oly+1,sNy+Oly
# Line 652  C       Ignore beta in the calculation o Line 824  C       Ignore beta in the calculation o
824            ENDDO            ENDDO
825           ENDDO           ENDDO
826          ENDDO          ENDDO
827            
828        ELSE        ELSE
829  C       Do not ignore beta  C       Do not ignore beta
830          DO k=1,Nr          DO k=1,Nr
# Line 679  C     This is necessary to find the eige Line 851  C     This is necessary to find the eige
851       &                *( N2(i,j,k)+N2(i,j-1,k) )       &                *( N2(i,j,k)+N2(i,j-1,k) )
852          ENDDO          ENDDO
853         ENDDO         ENDDO
 C     This fudge is necessary to avoid division by zero in gmredi_calc_eigs.  
 C     It does not affect the end result since it's in the overlap region.  
        j=1-Oly  
        DO i=1-Olx,sNx+Olx  
         N2W(i,j,k) = GM_K3D_minN2  
         N2S(i,j,k) = GM_K3D_minN2  
        ENDDO  
        i=1-Olx  
        DO j=1-Oly,sNy+Oly  
         N2W(i,j,k) = GM_K3D_minN2  
         N2S(i,j,k) = GM_K3D_minN2  
        ENDDO  
854        ENDDO        ENDDO
855    
856        IF(GM_K3D_likeGM) THEN  C     If GM_K3D_use_constK=.TRUE., the diffusivity for the eddy transport is
857    C     set to a constant equal to GM_K3D_constK.
858    C     If the diffusivity is constant the method here is the same as GM.
859    C     If GM_K3D_constRedi=.TRUE. K3D will be set equal to GM_K3D_constK later.
860          IF(GM_K3D_use_constK) THEN
861          DO k=1,Nr          DO k=1,Nr
862           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
863            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
# Line 701  C     It does not affect the end result Line 865  C     It does not affect the end result
865            ENDDO            ENDDO
866           ENDDO           ENDDO
867          ENDDO          ENDDO
       ELSE  
         DO k=1,Nr  
          DO j=1-Oly,sNy+Oly  
           DO i=1-Olx,sNx+Olx  
            KPV(i,j,k) = K3D(i,j,k,bi,bj)  
           ENDDO  
          ENDDO  
         ENDDO  
868        ENDIF        ENDIF
869    
870        IF (.NOT. GM_K3D_smooth) THEN        IF (.NOT. GM_K3D_smooth) THEN
871  C       Do not expand K grad(q) => no smoothing  C       Do not expand K grad(q) => no smoothing
872  C       May only be done with a constant K, otherwise the  C       May only be done with a constant K, otherwise the
873  C       integral constraint is violated.  C       integral constraint is violated.
874          DO k=1,Nr          DO k=1,Nr
875           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
# Line 739  C       Calculate the eigenvectors at th Line 895  C       Calculate the eigenvectors at th
895       I         rLowW(:,:,bi,bj),GM_K3D_NModes,.FALSE.,       I         rLowW(:,:,bi,bj),GM_K3D_NModes,.FALSE.,
896       O         dummy,modesW(:,:,:,:,bi,bj))       O         dummy,modesW(:,:,:,:,bi,bj))
897          ENDIF          ENDIF
898            
899  C       Calculate Xi_m at the west face of a cell  C       Calculate Xi_m at the west face of a cell
900          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
901           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 752  C       Calculate Xi_m at the west face Line 908  C       Calculate Xi_m at the west face
908           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
909            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
910             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
911                XimX(m,i,j) = XimX(m,i,j)              Kdqdx(i,j,k) = KPV(i,j,k)*gradqx(i,j,k)
912       &             - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)              XimX(m,i,j) = XimX(m,i,j)
913       &             *KPV(i,j,k)*gradqx(i,j,k)*modesW(m,i,j,k,bi,bj)       &           - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)
914         &           *Kdqdx(i,j,k)*modesW(m,i,j,k,bi,bj)
915             ENDDO             ENDDO
916            ENDDO            ENDDO
917           ENDDO           ENDDO
918          ENDDO          ENDDO
919              
920  C     Calculate Xi in the X direction at the west face  C     Calculate Xi in the X direction at the west face
921          DO k=1,Nr          DO k=1,Nr
922           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
# Line 772  C     Calculate Xi in the X direction at Line 929  C     Calculate Xi in the X direction at
929           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
930            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
931             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
932              Xix(i,j,k) = Xix(i,j,k)              Xix(i,j,k) = Xix(i,j,k)
933       &           + maskW(i,j,k,bi,bj)*XimX(m,i,j)*modesW(m,i,j,k,bi,bj)       &           + maskW(i,j,k,bi,bj)*XimX(m,i,j)*modesW(m,i,j,k,bi,bj)
934             ENDDO             ENDDO
935            ENDDO            ENDDO
# Line 802  C     Calculate the eigenvectors at the Line 959  C     Calculate the eigenvectors at the
959           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
960            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
961             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
962                Kdqdy(i,j,k) = KPV(i,j,k)*gradqy(i,j,k)
963              XimY(m,i,j) = XimY(m,i,j)              XimY(m,i,j) = XimY(m,i,j)
964       &           - drF(k)*hfacS(i,j,k,bi,bj)       &           - drF(k)*hfacS(i,j,k,bi,bj)
965       &           *KPV(i,j,k)*gradqy(i,j,k)*modesS(m,i,j,k,bi,bj)       &           *Kdqdy(i,j,k)*modesS(m,i,j,k,bi,bj)
966             ENDDO             ENDDO
967            ENDDO            ENDDO
968           ENDDO           ENDDO
969          ENDDO          ENDDO
970            
971  C     Calculate Xi for Y direction at the south face  C     Calculate Xi for Y direction at the south face
972          DO k=1,Nr          DO k=1,Nr
973           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
# Line 822  C     Calculate Xi for Y direction at th Line 980  C     Calculate Xi for Y direction at th
980           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
981            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
982             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
983              Xiy(i,j,k) = Xiy(i,j,k)              Xiy(i,j,k) = Xiy(i,j,k)
984       &           + maskS(i,j,k,bi,bj)*XimY(m,i,j)*modesS(m,i,j,k,bi,bj)       &           + maskS(i,j,k,bi,bj)*XimY(m,i,j)*modesS(m,i,j,k,bi,bj)
985             ENDDO             ENDDO
986            ENDDO            ENDDO
987           ENDDO           ENDDO
988          ENDDO          ENDDO
989    
990  C     ENDIF GM_K3D_likeGM  C     ENDIF (.NOT. GM_K3D_smooth)
991        ENDIF        ENDIF
992    
993    C     Calculate the renormalisation factor
994          DO j=1-Oly,sNy+Oly
995           DO i=1-Olx,sNx+Olx
996            uInt(i,j)=zeroRL
997            vInt(i,j)=zeroRL
998            KdqdyInt(i,j)=zeroRL
999            KdqdxInt(i,j)=zeroRL
1000            uKdqdyInt(i,j)=zeroRL
1001            vKdqdxInt(i,j)=zeroRL
1002            uXiyInt(i,j)=zeroRL
1003            vXixInt(i,j)=zeroRL
1004            Renorm(i,j)=oneRL
1005            RenormU(i,j)=oneRL
1006            RenormV(i,j)=oneRL
1007           ENDDO
1008          ENDDO
1009          DO k=1,Nr
1010           DO j=1-Oly,sNy+Oly-1
1011            DO i=1-Olx,sNx+Olx-1
1012             centreX = op5*(uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))
1013             centreY = op5*(Kdqdy(i,j,k)     +Kdqdy(i,j+1,k)     )
1014    C        For the numerator
1015             uInt(i,j) = uInt(i,j)
1016         &        + centreX*hfacC(i,j,k,bi,bj)*drF(k)
1017             KdqdyInt(i,j) = KdqdyInt(i,j)
1018         &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)
1019             uKdqdyInt(i,j) = uKdqdyInt(i,j)
1020         &        + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
1021    C        For the denominator
1022             centreY = op5*(Xiy(i,j,k) + Xiy(i,j+1,k))
1023             uXiyInt(i,j) = uXiyInt(i,j)
1024         &        + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
1025    
1026             centreX = op5*(Kdqdx(i,j,k)     +Kdqdx(i+1,j,k))
1027             centreY = op5*(vVel(i,j,k,bi,bj)+vVel(i,j+1,k,bi,bj) )
1028    C        For the numerator
1029             vInt(i,j) = vInt(i,j)
1030         &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)
1031             KdqdxInt(i,j) = KdqdxInt(i,j)
1032         &        + CentreX*hfacC(i,j,k,bi,bj)*drF(k)
1033             vKdqdxInt(i,j) = vKdqdxInt(i,j)
1034         &        + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
1035    C        For the denominator
1036             centreX = op5*(Xix(i,j,k) + Xix(i+1,j,k))
1037             vXixInt(i,j) = vXixInt(i,j)
1038         &        + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
1039    
1040            ENDDO
1041           ENDDO
1042          ENDDO
1043    
1044          DO j=1-Oly,sNy+Oly-1
1045           DO i=1-Olx,sNx+Olx-1
1046            IF (kLowC(i,j,bi,bj).GT.0) THEN
1047              numerator =
1048         &         (uKdqdyInt(i,j)-uInt(i,j)*KdqdyInt(i,j)/R_low(i,j,bi,bj))
1049         &        -(vKdqdxInt(i,j)-vInt(i,j)*KdqdxInt(i,j)/R_low(i,j,bi,bj))
1050              denominator = uXiyInt(i,j) - vXixInt(i,j)
1051    C         We can have troubles with floating point exceptions if the denominator
1052    C         of the renormalisation if the ocean is resting (e.g. intial conditions).
1053    C         So we make the renormalisation factor one if the denominator is very small
1054    C         The renormalisation factor is supposed to correct the error in the extraction of
1055    C         potential energy associated with the truncation of the expansion. Thus, we
1056    C         enforce a minimum value for the renormalisation factor.
1057    C         We also enforce a maximum renormalisation factor.
1058              IF (denominator.GT.small) THEN
1059                Renorm(i,j) = ABS(numerator/denominator)
1060                Renorm(i,j) = MAX(Renorm(i,j),GM_K3D_minRenorm)
1061                Renorm(i,j) = MIN(Renorm(i,j),GM_K3D_maxRenorm)
1062              ENDIF
1063            ENDIF
1064           ENDDO
1065          ENDDO
1066    C     Now put it back on to the velocity grids
1067          DO j=1-Oly+1,sNy+Oly-1
1068           DO i=1-Olx+1,sNx+Olx-1
1069            RenormU(i,j) = op5*(Renorm(i-1,j)+Renorm(i,j))
1070            RenormV(i,j) = op5*(Renorm(i,j-1)+Renorm(i,j))
1071           ENDDO
1072          ENDDO
1073    
1074  C     Calculate the eddy induced velocity in the X direction at the west face  C     Calculate the eddy induced velocity in the X direction at the west face
1075        DO k=1,Nr        DO k=1,Nr
1076         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
1077          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
1078           ustar(i,j,k) = -Xix(i,j,k)/coriU(i,j)           ustar(i,j,k) = -RenormU(i,j)*Xix(i,j,k)/coriU(i,j)
1079          ENDDO          ENDDO
1080         ENDDO         ENDDO
1081        ENDDO              ENDDO
1082    
1083  C     Calculate the eddy induced velocity in the Y direction at the south face  C     Calculate the eddy induced velocity in the Y direction at the south face
1084        DO k=1,Nr        DO k=1,Nr
1085         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
1086          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
1087           vstar(i,j,k) = -Xiy(i,j,k)/coriV(i,j)           vstar(i,j,k) = -RenormV(i,j)*Xiy(i,j,k)/coriV(i,j)
1088          ENDDO          ENDDO
1089         ENDDO         ENDDO
1090        ENDDO              ENDDO
1091    
1092  C     ======================================  C     ======================================
1093  C     Calculate the eddy induced overturning streamfunction  C     Calculate the eddy induced overturning streamfunction
# Line 869  C     ================================== Line 1107  C     ==================================
1107          ENDDO          ENDDO
1108         ENDDO         ENDDO
1109        ENDDO        ENDDO
1110        
1111  #else  #else
1112    
1113        IF (GM_AdvForm) THEN        IF (GM_AdvForm) THEN
# Line 897  C     ================================== Line 1135  C     ==================================
1135  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
1136  C     Diagnostics  C     Diagnostics
1137        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
1138          CALL DIAGNOSTICS_FILL(K3D,    'GM_K3D  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(K3D,    'GM_K3D  ',0,Nr,1,bi,bj,myThid)
1139          CALL DIAGNOSTICS_FILL(urms,   'GM_URMS ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(KPV,    'GM_KPV  ',0,Nr,2,bi,bj,myThid)
1140          CALL DIAGNOSTICS_FILL(Rdef,   'GM_RDEF ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(urms,   'GM_URMS ',0,Nr,2,bi,bj,myThid)
1141          CALL DIAGNOSTICS_FILL(RRhines,'GM_LRHNS',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Rdef,   'GM_RDEF ',0, 1,1,bi,bj,myThid)
1142          CALL DIAGNOSTICS_FILL(Rmix,   'GM_RMIX ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Rurms,  'GM_RURMS',0, 1,2,bi,bj,myThid)
1143          CALL DIAGNOSTICS_FILL(supp,   'GM_SUPP ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(RRhines,'GM_RRHNS',0, 1,2,bi,bj,myThid)
1144          CALL DIAGNOSTICS_FILL(Xix,    'GM_Xix  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Rmix,   'GM_RMIX ',0, 1,2,bi,bj,myThid)
1145          CALL DIAGNOSTICS_FILL(Xiy,    'GM_Xiy  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(supp,   'GM_SUPP ',0,Nr,2,bi,bj,myThid)
1146          CALL DIAGNOSTICS_FILL(cDopp,  'GM_C    ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Xix,    'GM_Xix  ',0,Nr,2,bi,bj,myThid)
1147          CALL DIAGNOSTICS_FILL(Ubaro,  'GM_UBARO',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Xiy,    'GM_Xiy  ',0,Nr,2,bi,bj,myThid)
1148          CALL DIAGNOSTICS_FILL(eady,   'GM_EADY ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(cDopp,  'GM_C    ',0, 1,2,bi,bj,myThid)
1149          CALL DIAGNOSTICS_FILL(SlopeX, 'GM_Sx   ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Ubaro,  'GM_UBARO',0, 1,2,bi,bj,myThid)
1150          CALL DIAGNOSTICS_FILL(SlopeY, 'GM_Sy   ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(eady,   'GM_EADY ',0, 1,2,bi,bj,myThid)
1151          CALL DIAGNOSTICS_FILL(tfluxX, 'GM_TFLXX',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(SlopeX, 'GM_Sx   ',0,Nr,2,bi,bj,myThid)
1152          CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(SlopeY, 'GM_Sy   ',0,Nr,2,bi,bj,myThid)
1153          CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(tfluxX, 'GM_TFLXX',0,Nr,2,bi,bj,myThid)
1154          CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,2,bi,bj,myThid)
1155          CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,2,bi,bj,myThid)
1156          CALL DIAGNOSTICS_FILL(ustar,  'GM_USTAR',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,2,bi,bj,myThid)
1157          CALL DIAGNOSTICS_FILL(vstar,  'GM_VSTAR',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Kdqdy,  'GM_Kdqdy',0,Nr,2,bi,bj,myThid)
1158          CALL DIAGNOSTICS_FILL(umc,    'GM_UMC  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Kdqdx,  'GM_Kdqdx',0,Nr,2,bi,bj,myThid)
1159          CALL DIAGNOSTICS_FILL(ubar,   'GM_UBAR ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,2,bi,bj,myThid)
1160          CALL DIAGNOSTICS_FILL(modesC(1,:,:,:,bi,bj),          CALL DIAGNOSTICS_FILL(ustar,  'GM_USTAR',0,Nr,2,bi,bj,myThid)
1161       &                                'GM_MODEC',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(vstar,  'GM_VSTAR',0,Nr,2,bi,bj,myThid)
1162            CALL DIAGNOSTICS_FILL(umc,    'GM_UMC  ',0,Nr,2,bi,bj,myThid)
1163            CALL DIAGNOSTICS_FILL(ubar,   'GM_UBAR ',0,Nr,2,bi,bj,myThid)
1164            CALL DIAGNOSTICS_FILL(modesC, 'GM_MODEC',0,Nr,1,bi,bj,myThid)
1165            CALL DIAGNOSTICS_FILL(M4loc,  'GM_M4   ',0,Nr,2,bi,bj,myThid)
1166            CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,2,bi,bj,myThid)
1167            CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,2,bi,bj,myThid)
1168            CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,2,bi,bj,myThid)
1169            CALL DIAGNOSTICS_FILL(Renorm, 'GM_RENRM',0, 1,2,bi,bj,myThid)
1170    #ifdef GM_K3D_PASSIVE
1171            CALL DIAGNOSTICS_FILL(psistar,'GM_PSTAR',0,Nr,2,bi,bj,myThid)
1172    #endif
1173        ENDIF        ENDIF
1174  #endif  #endif
1175    
1176    C     For the Redi diffusivity, we set K3D to a constant if
1177    C     GM_K3D_constRedi=.TRUE.
1178          IF (GM_K3D_constRedi) THEN
1179            DO k=1,Nr
1180             DO j=1-Oly,sNy+Oly
1181              DO i=1-Olx,sNx+Olx
1182               K3D(i,j,k,bi,bj) = GM_K3D_constK
1183              ENDDO
1184             ENDDO
1185            ENDDO
1186          ENDIF
1187    
1188    #ifdef ALLOW_DIAGNOSTICS
1189          IF ( useDiagnostics )
1190         &     CALL DIAGNOSTICS_FILL(K3D, 'GM_K3D_T',0,Nr,1,bi,bj,myThid)
1191    #endif
1192    
1193  #endif /* GM_K3D */  #endif /* GM_K3D */
1194        RETURN        RETURN
1195        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22