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

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

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

Revision 1.14 - (hide annotations) (download)
Wed Jan 1 23:20:48 2014 UTC (10 years, 5 months ago) by m_bates
Branch: MAIN
CVS Tags: checkpoint64s, checkpoint64u, checkpoint64t
Changes since 1.13: +99 -78 lines 
K3d: Added option to keep diffusivity constant in the mixed layer. Also added some diagnostics.
1 m_bates 1.14 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_k3d.F,v 1.13 2013/10/21 18:46:05 m_bates Exp $
2 m_bates 1.1 C $Name: $
3 m_bates 1.9 #include "CPP_OPTIONS.h"
4 m_bates 1.1 #include "GMREDI_OPTIONS.h"
5    
6     C !ROUTINE: GMREDI_K3D
7     C !INTERFACE:
8     SUBROUTINE GMREDI_K3D(
9     I iMin, iMax, jMin, jMax,
10     I sigmaX, sigmaY, sigmaR,
11 m_bates 1.4 I bi, bj, myTime, myThid )
12 m_bates 1.1
13     C !DESCRIPTION: \bv
14     C *==========================================================*
15     C | SUBROUTINE GMREDI_K3D
16     C | o Calculates the 3D diffusivity as per Bates et al. (2013)
17     C *==========================================================*
18     C \ev
19    
20     IMPLICIT NONE
21    
22     C == Global variables ==
23     #include "SIZE.h"
24     #include "GRID.h"
25     #include "DYNVARS.h"
26     #include "EEPARAMS.h"
27     #include "PARAMS.h"
28     #include "GMREDI.h"
29    
30     C !INPUT/OUTPUT PARAMETERS:
31     C == Routine arguments ==
32     C bi, bj :: tile indices
33     C myThid :: My Thread Id. number
34    
35     INTEGER iMin,iMax,jMin,jMax
36     _RL sigmaX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
37     _RL sigmaY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
38     _RL sigmaR(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
39     INTEGER bi, bj
40 m_bates 1.4 _RL myTime
41 m_bates 1.1 INTEGER myThid
42    
43     #ifdef GM_K3D
44    
45 m_bates 1.4 C === Functions ====
46     LOGICAL DIFFERENT_MULTIPLE
47     EXTERNAL DIFFERENT_MULTIPLE
48    
49 m_bates 1.1 C !LOCAL VARIABLES:
50     C == Local variables ==
51 m_bates 1.4 INTEGER i,j,k,kk,m
52    
53     C update_modes :: Whether to update the eigenmodes
54     LOGICAL update_modes
55    
56     C surfk :: index of the depth of the surface layer
57     C kLow_C :: Local version of the index of deepest wet grid cell on tracer grid
58     C kLow_U :: Local version of the index of deepest wet grid cell on U grid
59     C kLow_V :: Local version of the index of deepest wet grid cell on V grid
60 m_bates 1.1 INTEGER surfk(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
61     INTEGER kLow_C(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
62     INTEGER kLow_U(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
63     INTEGER kLow_V(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
64 m_bates 1.4
65     C N2loc :: local N**2
66     C slope :: local slope
67     C Req :: local equatorial deformation radius (m)
68     C deltaH :: local thickness of Eady integration (m)
69     C g_reciprho_sq :: (gravity*recip_rhoConst)**2
70     C M4loc :: local M**4
71     C maxDRhoDz :: maximum value of d(rho)/dz (derived from GM_K3D_minN2)
72     C sigx :: local d(rho)/dx
73     C sigy :: local d(rho)/dy
74     C sigz :: local d(rho)/dz
75     C hsurf :: local surface layer depth
76     C small :: a small number (to avoid floating point exceptions)
77 m_bates 1.10 _RL N2loc(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
78 m_bates 1.1 _RL slope
79 m_bates 1.10 _RL slopeC(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
80 m_bates 1.4 _RL Req
81 m_bates 1.10 _RL deltaH(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82 m_bates 1.1 _RL g_reciprho_sq
83 m_bates 1.10 _RL M4loc(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
84 m_bates 1.13 _RL M4onN2(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
85 m_bates 1.1 _RL maxDRhoDz
86     _RL sigx, sigy, sigz
87     _RL hsurf
88     _RL small
89 m_bates 1.4
90     C dfdy :: gradient of the Coriolis paramter, df/dy, 1/(m*s)
91     C dfdx :: gradient of the Coriolis paramter, df/dx, 1/(m*s)
92 m_bates 1.6 C gradf :: gradient of the Coriolis paramter at a cell centre, 1/(m*s)
93 m_bates 1.8 C Rurms :: a local mixing length used in calculation of urms (m)
94 m_bates 1.4 C RRhines :: The Rhines scale (m)
95     C Rmix :: Mixing length
96     C N2 :: Square of the buoyancy frequency (1/s**2)
97     C N2W :: Square of the buoyancy frequency (1/s**2) averaged to west of grid cell
98     C N2S :: Square of the buoyancy frequency (1/s**2) averaged to south of grid cell
99     C N :: Buoyancy frequency, SQRT(N2)
100     C BVint :: The vertical integral of N (m/s)
101     C ubar :: Zonal velocity on a tracer point (m/s)
102     C vbar :: Meridional velocity on a tracer point (m/s)
103     C Ubaro :: Barotropic velocity (m/s)
104     _RL dfdy( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
105     _RL dfdx( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
106     _RL gradf( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
107     _RL dummy( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
108 m_bates 1.8 _RL Rurms( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
109 m_bates 1.4 _RL RRhines(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
110     _RL Rmix( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
111     _RL N2( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
112     _RL N2W( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
113     _RL N2S( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
114     _RL N( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
115     _RL BVint( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
116     _RL Ubaro( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
117     _RL ubar( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
118     _RL vbar( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
119    
120     C Rmid :: Rossby radius (m)
121     C KPV :: Diffusivity (m**2/s)
122 m_bates 1.13 C Kdqdx :: diffusivity multiplied by zonal PV gradient
123     C Kdqdy :: diffusivity multiplied by meridional PV gradient
124 m_bates 1.4 C SlopeX :: isopycnal slope in x direction
125     C SlopeY :: isopycnal slope in y direction
126     C dSigmaDx :: sigmaX averaged onto tracer grid
127     C dSigmaDy :: sigmaY averaged onto tracer grid
128     C tfluxX :: thickness flux in x direction
129     C tfluxY :: thickness flux in y direction
130     C fCoriU :: Coriolis parameter averaged to U points
131     C fCoriV :: Coriolis parameter averaged to V points
132     C cori :: Coriolis parameter forced to be finite near the equator
133     C coriU :: As for cori, but, at U point
134     C coriV :: As for cori, but, at V point
135     C surfkz :: Depth of surface layer (in r units)
136     _RL Rmid(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
137 m_bates 1.3 _RL KPV(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
138 m_bates 1.13 _RL Kdqdy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
139     _RL Kdqdx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
140 m_bates 1.1 _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
141     _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
142     _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
143     _RL dSigmaDy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
144     _RL tfluxX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
145     _RL tfluxY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
146     _RL cori(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
147     _RL coriU(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
148     _RL coriV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
149     _RL fCoriU(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
150     _RL fCoriV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
151     _RL surfkz(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
152 m_bates 1.4
153     C gradqx :: Potential vorticity gradient in x direction
154     C gradqy :: Potential vorticity gradient in y direction
155     C XimX :: Vertical integral of phi_m*K*gradqx
156     C XimY :: Vertical integral of phi_m*K*gradqy
157     C cDopp :: Quasi-Doppler shifted long Rossby wave speed (m/s)
158     C umc :: ubar-c (m/s)
159     C eady :: Eady growth rate (1/s)
160     C urms :: the rms eddy velocity (m/s)
161     C supp :: The suppression factor
162     C ustar :: The eddy induced velocity in the x direction
163     C vstar :: The eddy induced velocity in the y direction
164     C Xix :: Xi in the x direction (m/s**2)
165     C Xiy :: Xi in the y direction (m/s**2)
166 m_bates 1.1 _RL gradqx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
167     _RL gradqy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
168     _RL XimX(GM_K3D_NModes,1-Olx:sNx+Olx,1-Oly:sNy+Oly)
169     _RL XimY(GM_K3D_NModes,1-Olx:sNx+Olx,1-Oly:sNy+Oly)
170 m_bates 1.4 _RL cDopp(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
171 m_bates 1.1 _RL umc( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
172     _RL eady( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
173     _RL urms( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
174     _RL supp( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
175     _RL ustar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
176     _RL vstar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
177     _RL Xix( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
178     _RL Xiy( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
179 m_bates 1.4 #ifdef GM_K3D_PASSIVE
180     C psistar :: eddy induced streamfunction in the y direction
181 m_bates 1.1 _RL psistar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
182 m_bates 1.4 #endif
183 m_bates 1.1
184    
185     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
186    
187     C ======================================
188     C Initialise some variables
189     C ======================================
190     small = TINY(zeroRL)
191 m_bates 1.4 update_modes=.FALSE.
192     IF ( DIFFERENT_MULTIPLE(GM_K3D_vecFreq,myTime,deltaTClock) )
193     & update_modes=.TRUE.
194    
195 m_bates 1.1 DO j=1-Oly,sNy+Oly
196     DO i=1-Olx,sNx+Olx
197     kLow_C(i,j) = kLowC(i,j,bi,bj)
198     ENDDO
199     ENDDO
200     DO j=1-Oly,sNy+Oly
201     DO i=1-Olx+1,sNx+Olx
202     kLow_U(i,j) = MIN( kLow_C(i,j), kLow_C(i-1,j) )
203     ENDDO
204     ENDDO
205     DO j=1-Oly+1,sNy+Oly
206     DO i=1-Olx,sNx+Olx
207     kLow_V(i,j) = MIN( kLow_C(i,j), kLow_C(i,j-1) )
208     ENDDO
209     ENDDO
210    
211     C Dummy values for the edges
212     C This avoids weirdness in gmredi_calc_eigs
213     i=1-Olx
214     DO j=1-Oly,sNy+Oly
215     kLow_U(i,j) = kLow_C(i,j)
216     ENDDO
217     j=1-Oly
218     DO i=1-Olx,sNx+Olx
219     kLow_V(i,j) = kLow_C(i,j)
220     ENDDO
221    
222     g_reciprho_sq = (gravity*recip_rhoConst)**2
223     C Gradient of Coriolis
224     DO j=1-Oly+1,sNy+Oly
225     DO i=1-Olx+1,sNx+Olx
226     dfdx(i,j) = ( fCori(i,j,bi,bj)-fCori(i-1,j,bi,bj) )
227     & *recip_dxC(i,j,bi,bj)
228     dfdy(i,j) = ( fCori(i,j,bi,bj)-fCori(i,j-1,bi,bj) )
229     & *recip_dyC(i,j,bi,bj)
230     ENDDO
231     ENDDO
232    
233     C Coriolis at C points enforcing a minimum value so
234     C that it is defined at the equator
235     DO j=1-Oly,sNy+Oly
236     DO i=1-Olx,sNx+Olx
237     cori(i,j) = SIGN( MAX( ABS(fCori(i,j,bi,bj)),GM_K3D_minCori ),
238     & fCori(i,j,bi,bj) )
239     ENDDO
240     ENDDO
241     C Coriolis at U and V points
242     DO j=1-Oly+1,sNy+Oly
243     DO i=1-Olx+1,sNx+Olx
244     C Limited so that the inverse is defined at the equator
245     coriU(i,j) = op5*( cori(i,j)+cori(i-1,j) )
246     coriU(i,j) = SIGN( MAX( ABS(coriU(i,j)),GM_K3D_minCori ),
247     & coriU(i,j) )
248    
249     coriV(i,j) = op5*( cori(i,j)+cori(i,j-1) )
250     coriV(i,j) = SIGN( MAX( ABS(coriV(i,j)),GM_K3D_minCori ),
251     & coriV(i,j) )
252    
253     C Not limited
254     fCoriU(i,j) = op5*( fCori(i,j,bi,bj)+fCori(i-1,j,bi,bj) )
255     fCoriV(i,j) = op5*( fCori(i,j,bi,bj)+fCori(i,j-1,bi,bj) )
256     ENDDO
257     ENDDO
258     DO j=1-Oly,sNy+Oly
259     DO i=1-Olx,sNx+Olx
260 m_bates 1.6 gradf(i,j) = recip_rSphere*fCoriCos(i,j,bi,bj)
261 m_bates 1.1 ENDDO
262     ENDDO
263    
264     C Zeroing some cumulative fields
265     DO j=1-Oly,sNy+Oly
266     DO i=1-Olx,sNx+Olx
267 m_bates 1.10 eady(i,j) = zeroRL
268     BVint(i,j) = zeroRL
269     Ubaro(i,j) = zeroRL
270     deltaH(i,j) = zeroRL
271     ENDDO
272     ENDDO
273     DO k=1,Nr
274     DO j=1-Oly,sNy+Oly
275     DO i=1-Olx,sNx+Olx
276     slopeC(i,j,k)=zeroRL
277     ENDDO
278 m_bates 1.1 ENDDO
279     ENDDO
280    
281     C Find the zonal velocity at the cell centre
282     C The logicals here are, in order: 1/ go from grid to north/east directions
283     C 2/ go from C to A grid and 3/ apply the mask
284 m_bates 1.9 #ifdef ALLOW_EDDYPSI
285     IF (GM_InMomAsStress) THEN
286     CALL rotate_uv2en_rl(uMean, vMean, ubar, vbar, .TRUE., .TRUE.,
287     & .TRUE.,Nr,mythid)
288     ELSE
289     #endif
290     CALL rotate_uv2en_rl(uVel, vVel, ubar, vbar, .TRUE., .TRUE.,
291     & .TRUE.,Nr,mythid)
292     #ifdef ALLOW_EDDYPSI
293     ENDIF
294     #endif
295 m_bates 1.1
296     C Square of the buoyancy frequency at the top of a grid cell
297     DO k=2,Nr
298     DO j=1-Oly,sNy+Oly
299     DO i=1-Olx,sNx+Olx
300     N2(i,j,k) = -gravity*recip_rhoConst*sigmaR(i,j,k)
301     ENDDO
302     ENDDO
303     ENDDO
304     C N2(k=1) is always zero
305     k=1
306     DO j=1-Oly,sNy+Oly
307     DO i=1-Olx,sNx+Olx
308     N2(i,j,k) = 0.0
309     N(i,j,k) = 0.0
310     ENDDO
311     ENDDO
312     C Enforce a minimum N2
313     DO k=2,Nr
314     DO j=1-Oly,sNy+Oly
315     DO i=1-Olx,sNx+Olx
316     IF (N2(i,j,k).LT.GM_K3D_minN2) N2(i,j,k)=GM_K3D_minN2
317     N(i,j,k) = SQRT(N2(i,j,k))
318     ENDDO
319     ENDDO
320     ENDDO
321     C Calculate the minimum drho/dz
322     maxDRhoDz = -rhoConst*GM_K3D_minN2/gravity
323    
324     C Calculate the barotropic velocity by vertically integrating
325     C and the dividing by the depth of the water column
326     C Note that Ubaro is on the U grid.
327     DO k=1,Nr
328     DO j=1-Oly,sNy+Oly
329     DO i=1-Olx,sNx+Olx
330     Ubaro(i,j) = Ubaro(i,j) +
331     & maskW(i,j,k,bi,bj)*drF(k)*hfacC(i,j,k,bi,bj)
332     & *ubar(i,j,k,bi,bj)
333     ENDDO
334     ENDDO
335     ENDDO
336     DO j=1-Oly,sNy+Oly
337     DO i=1-Olx,sNx+Olx
338     IF (kLow_C(i,j).GT.0) THEN
339     C The minus sign is because r_Low<0
340     Ubaro(i,j) = -Ubaro(i,j)/r_Low(i,j,bi,bj)
341     ENDIF
342     ENDDO
343     ENDDO
344    
345     C Integrate the buoyancy frequency vertically using the trapezoidal method.
346     DO k=1,Nr
347     DO j=1-Oly,sNy+Oly
348     DO i=1-Olx,sNx+Olx
349     IF (k.LT.kLow_C(i,j)) THEN
350     BVint(i,j) = BVint(i,j) + hFacC(i,j,k,bi,bj)*drF(k)
351     & *(N(i,j,k)+N(i,j,k+1))
352     ELSEIF (k.EQ.kLow_C(i,j)) THEN
353     C Assume that N(z=-H)=0
354     BVint(i,j) = BVint(i,j) + hFacC(i,j,k,bi,bj)*drF(k)*N(i,j,k)
355     ENDIF
356     ENDDO
357     ENDDO
358     ENDDO
359     DO j=1-Oly,sNy+Oly
360     DO i=1-Olx,sNx+Olx
361     BVint(i,j) = op5*BVint(i,j)
362     ENDDO
363     ENDDO
364    
365     C Calculate the eigenvalues and eigenvectors
366 m_bates 1.4 IF (update_modes) THEN
367     CALL GMREDI_CALC_EIGS(
368     I iMin,iMax,jMin,jMax,bi,bj,N2,myThid,
369     I kLow_C, maskC(:,:,:,bi,bj),
370     I hfacC(:,:,:,bi,bj), recip_hfacC(:,:,:,bi,bj),
371     I R_Low(:,:,bi,bj), 1, .TRUE.,
372     O Rmid, modesC(:,:,:,:,bi,bj))
373    
374     C Calculate the Rossby Radius
375     DO j=1-Oly+1,sNy+Oly
376     DO i=1-Olx+1,sNx+Olx
377     Req = SQRT(BVint(i,j)/(2*pi*gradf(i,j)))
378     Rdef(i,j,bi,bj) = MIN(Rmid(i,j),Req)
379     ENDDO
380     ENDDO
381     ENDIF
382 m_bates 1.1
383     C Average dsigma/dx and dsigma/dy onto the centre points
384    
385     DO k=1,Nr
386     DO j=1-Oly,sNy+Oly-1
387     DO i=1-Olx,sNx+Olx-1
388     dSigmaDx(i,j,k) = op5*(sigmaX(i,j,k)+sigmaX(i+1,j,k))
389     dSigmaDy(i,j,k) = op5*(sigmaY(i,j,k)+sigmaY(i,j+1,k))
390     ENDDO
391     ENDDO
392     ENDDO
393    
394     C ===============================
395     C Calculate the Eady growth rate
396     C ===============================
397     DO k=1,Nr
398    
399 m_bates 1.10 DO j=1-Oly,sNy+Oly-1
400     DO i=1-Olx,sNx+Olx-1
401     M4loc(i,j,k) = g_reciprho_sq*( dSigmaDx(i,j,k)**2
402     & +dSigmaDy(i,j,k)**2 )
403     IF (k.NE.kLow_C(i,j)) THEN
404     N2loc(i,j,k) = op5*(N2(i,j,k)+N2(i,j,k+1))
405     ELSE
406     N2loc(i,j,k) = op5*N2(i,j,k)
407     ENDIF
408     ENDDO
409     ENDDO
410 m_bates 1.1 C The bottom of the grid cell is shallower than the top
411     C integration level, so, advance the depth.
412 m_bates 1.10 IF (-rF(k+1) .LE. GM_K3D_EadyMinDepth) CYCLE
413 m_bates 1.1
414 jmc 1.7 C Do not bother going any deeper since the top of the
415 m_bates 1.1 C cell is deeper than the bottom integration level
416     IF (-rF(k).GE.GM_K3D_EadyMaxDepth) EXIT
417    
418     C We are in the integration depth range
419     DO j=1-Oly,sNy+Oly-1
420     DO i=1-Olx,sNx+Olx-1
421 m_bates 1.10 IF ( (kLow_C(i,j).GE.k) .AND.
422     & (-hMixLayer(i,j,bi,bj).LE.-rC(k)) ) THEN
423 m_bates 1.1
424 m_bates 1.12 slopeC(i,j,k) = SQRT(M4loc(i,j,k))/N2loc(i,j,k)
425 m_bates 1.1 C Limit the slope. Note, this is not all the Eady calculations.
426 m_bates 1.13 IF (slopeC(i,j,k).LE.GM_maxSlope) THEN
427     M4onN2(i,j,k) = M4loc(i,j,k)/N2loc(i,j,k)
428 m_bates 1.1 ELSE
429 m_bates 1.13 slopeC(i,j,k) = GM_maxslope
430     M4onN2(i,j,k) = SQRT(M4loc(i,j,k))*GM_maxslope
431 m_bates 1.1 ENDIF
432 m_bates 1.13 eady(i,j) = eady(i,j)
433     & + hfacC(i,j,k,bi,bj)*drF(k)*M4onN2(i,j,k)
434 m_bates 1.10 deltaH(i,j) = deltaH(i,j) + drF(k)
435 m_bates 1.1 ENDIF
436     ENDDO
437     ENDDO
438     ENDDO
439    
440     DO j=1-Oly,sNy+Oly
441     DO i=1-Olx,sNx+Olx
442 m_bates 1.10 C If the minimum depth for the integration is deeper than the ocean
443     C bottom OR the mixed layer is deeper than the maximum depth of
444     C integration, we set the Eady growth rate to something small
445     C to avoid floating point exceptions.
446     C Later, these areas will be given a small diffusivity.
447     IF (deltaH(i,j).EQ.zeroRL) THEN
448 m_bates 1.1 eady(i,j) = small
449    
450 m_bates 1.10 C Otherwise, divide over the integration and take the square root
451     C to actually find the Eady growth rate.
452 m_bates 1.1 ELSE
453 m_bates 1.10 eady(i,j) = SQRT(eady(i,j)/deltaH(i,j))
454 m_bates 1.1
455     ENDIF
456    
457     ENDDO
458     ENDDO
459    
460     C ======================================
461     C Calculate the diffusivity
462     C ======================================
463     DO j=1-Oly+1,sNy+Oly
464     DO i=1-Olx+1,sNx+Olx-1
465     C Calculate the Visbeck velocity
466 m_bates 1.13 Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_K3D_Rmax)
467 m_bates 1.8 urms(i,j,1) = GM_K3D_Lambda*eady(i,j)*Rurms(i,j)
468 m_bates 1.1 C Set the bottom urms to zero
469     k=kLow_C(i,j)
470     IF (k.GT.0) urms(i,j,k) = 0.0
471    
472     C Calculate the Rhines scale
473     RRhines(i,j) = SQRT(urms(i,j,1)/gradf(i,j))
474    
475     C Calculate the estimated length scale
476 m_bates 1.4 Rmix(i,j) = MIN(Rdef(i,j,bi,bj), RRhines(i,j))
477 m_bates 1.13 Rmix(i,j) = MAX(Rmix(i,j),GM_K3D_Rmin)
478 m_bates 1.1
479     C Calculate the Doppler shifted long Rossby wave speed
480     C Ubaro is on the U grid so we must average onto the M grid.
481     cDopp(i,j) = op5*( Ubaro(i,j)+Ubaro(i+1,j) )
482 m_bates 1.4 & - gradf(i,j)*Rdef(i,j,bi,bj)*Rdef(i,j,bi,bj)
483 m_bates 1.1 C Limit the wave speed to the namelist variable GM_K3D_maxC
484     IF (ABS(cDopp(i,j)).GT.GM_K3D_maxC) THEN
485     cDopp(i,j) = MAX(GM_K3D_maxC, cDopp(i,j))
486     ENDIF
487    
488     ENDDO
489     ENDDO
490    
491     C Project the surface urms to the subsurface using the first baroclinic mode
492 m_bates 1.4 CALL GMREDI_CALC_URMS(
493     I iMin,iMax,jMin,jMax,bi,bj,N2,myThid,
494     U urms)
495 m_bates 1.1
496     C Calculate the diffusivity (on the mass grid)
497     DO k=1,Nr
498     DO j=1-Oly,sNy+Oly
499     DO i=1-Olx,sNx+Olx
500     IF (k.LE.kLow_C(i,j)) THEN
501 m_bates 1.10 IF (deltaH(i,j).EQ.zeroRL) THEN
502 m_bates 1.1 K3D(i,j,k,bi,bj) = GM_K3D_smallK
503     ELSE
504     IF (urms(i,j,k).EQ.0.0) THEN
505     K3D(i,j,k,bi,bj) = GM_K3D_smallK
506     ELSE
507 m_bates 1.13 umc(i,j,k) =ubar(i,j,k,bi,bj) - cDopp(i,j)
508     supp(i,j,k)=1/(1+GM_K3D_b1*umc(i,j,k)**2/urms(i,j,1)**2)
509     C 2*Rmix gives the diameter
510 m_bates 1.1 K3D(i,j,k,bi,bj) = GM_K3D_gamma*urms(i,j,k)
511 m_bates 1.13 & *2*Rmix(i,j)*supp(i,j,k)
512 m_bates 1.1 ENDIF
513    
514     C Enforce lower and upper bounds on the diffusivity
515     IF (K3D(i,j,k,bi,bj).LT.GM_K3D_smallK)
516     & K3D(i,j,k,bi,bj) = GM_K3D_smallK
517     IF (K3D(i,j,k,bi,bj).GT.GM_maxK3D)
518     & K3D(i,j,k,bi,bj) = GM_maxK3D
519     ENDIF
520     ENDIF
521     ENDDO
522     ENDDO
523     ENDDO
524    
525     C ======================================
526     C Find the PV gradient
527     C ======================================
528 m_bates 1.3 C Calculate the surface layer thickness.
529     C Use hMixLayer (calculated in model/src/calc_oce_mxlayer)
530     C for the mixed layer depth.
531 m_bates 1.1
532 m_bates 1.3 C Enforce a minimum surface layer depth
533 m_bates 1.1 DO j=1-Oly,sNy+Oly
534     DO i=1-Olx,sNx+Olx
535 m_bates 1.3 surfkz(i,j) = MIN(-GM_K3D_surfMinDepth,-hMixLayer(i,j,bi,bj))
536     surfkz(i,j) = MAX(surfkz(i,j),R_low(i,j,bi,bj))
537     IF(maskC(i,j,1,bi,bj).EQ.0.0) surfkz(i,j)=0.0
538     surfk(i,j) = 0
539 m_bates 1.1 ENDDO
540     ENDDO
541 m_bates 1.4 DO k=1,Nr
542 m_bates 1.1 DO j=1-Oly,sNy+Oly
543     DO i=1-Olx,sNx+Olx
544 m_bates 1.3 IF (rF(k).GT.surfkz(i,j) .AND. surfkz(i,j).GE.rF(k+1))
545     & surfk(i,j) = k
546 m_bates 1.1 ENDDO
547     ENDDO
548     ENDDO
549 m_bates 1.3
550 m_bates 1.1 C Calculate the isopycnal slope
551     DO j=1-Oly,sNy+Oly-1
552     DO i=1-Olx,sNx+Olx-1
553     SlopeX(i,j,1) = zeroRL
554     SlopeY(i,j,1) = zeroRL
555     ENDDO
556     ENDDO
557     DO k=2,Nr
558     DO j=1-Oly+1,sNy+Oly
559     DO i=1-Olx+1,sNx+Olx
560     IF(surfk(i,j).GE.kLowC(i,j,bi,bj)) THEN
561     C If the surface layer is thinner than the water column
562     C the set the slope to zero to avoid problems.
563     SlopeX(i,j,k) = zeroRL
564     SlopeY(i,j,k) = zeroRL
565    
566     ELSE
567     C Calculate the zonal slope at the western cell face (U grid)
568 m_bates 1.4 sigz = MIN( op5*(sigmaR(i,j,k)+sigmaR(i-1,j,k)), maxDRhoDz )
569 m_bates 1.1 sigx = op5*( sigmaX(i,j,k)+sigmaX(i,j,k-1) )
570     slope = sigx/sigz
571     IF(ABS(slope).GT.GM_maxSlope)
572     & slope = SIGN(GM_maxSlope,slope)
573     SlopeX(i,j,k)=-maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)*slope
574    
575     C Calculate the meridional slope at the southern cell face (V grid)
576 m_bates 1.4 sigz = MIN( op5*(sigmaR(i,j,k)+sigmaR(i,j-1,k)), maxDRhoDz )
577 m_bates 1.1 sigy = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )
578     slope = sigy/sigz
579     IF(ABS(slope).GT.GM_maxSlope)
580     & slope = SIGN(GM_maxSlope,slope)
581     SlopeY(i,j,k)=-maskS(i,j,k-1,bi,bj)*maskS(i,j,k,bi,bj)*slope
582     ENDIF
583     ENDDO
584     ENDDO
585     ENDDO
586    
587 m_bates 1.14 C Calculate the thickness flux and diffusivity. These may be altered later
588     C depending on namelist options.
589 m_bates 1.1 C Enforce a zero slope bottom boundary condition for the bottom most cells (k=Nr)
590     k=Nr
591     DO j=1-Oly,sNy+Oly
592     DO i=1-Olx,sNx+Olx
593     C Zonal thickness flux at the western cell face
594     tfluxX(i,j,k) = -fCoriU(i,j)*SlopeX(i,j,k)
595     & *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
596     C Meridional thickness flux at the southern cell face
597     tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)
598     & *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
599 m_bates 1.14
600     C Use the interior diffusivity. Note: if we are using a
601     C constant diffusivity KPV is overwritten later
602     KPV(i,j,k) = K3D(i,j,k,bi,bj)
603    
604 m_bates 1.1 ENDDO
605     ENDDO
606    
607 m_bates 1.14 C Calculate the thickness flux and diffusivity and for other cells (k<Nr)
608 m_bates 1.1 DO k=Nr-1,1,-1
609 m_bates 1.14 DO j=1-Oly,sNy+Oly
610     DO i=1-Olx,sNx+Olx
611     C Zonal thickness flux at the western cell face
612     tfluxX(i,j,k)=-fCoriU(i,j)*(SlopeX(i,j,k)-SlopeX(i,j,k+1))
613     & *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
614     & *maskW(i,j,k,bi,bj)
615    
616     C Meridional thickness flux at the southern cell face
617     tfluxY(i,j,k)=-fCoriV(i,j)*(SlopeY(i,j,k)-SlopeY(i,j,k+1))
618     & *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
619     & *maskS(i,j,k,bi,bj)
620    
621     C Use the interior diffusivity. Note: if we are using a
622     C constant diffusivity KPV is overwritten later
623     KPV(i,j,k) = K3D(i,j,k,bi,bj)
624     ENDDO
625     ENDDO
626     ENDDO
627    
628    
629     C Special treatment for the surface layer (if necessary), which overwrites
630     C values in the previous loops.
631     IF (GM_K3D_ThickSheet .OR. GM_K3D_surfK) THEN
632     DO k=Nr-1,1,-1
633     DO j=1-Oly,sNy+Oly
634     DO i=1-Olx,sNx+Olx
635     IF(k.LE.surfk(i,j)) THEN
636     C We are in the surface layer. Change the thickness flux
637     C and diffusivity as necessary.
638    
639     IF (GM_K3D_ThickSheet) THEN
640     C We are in the surface layer, so set the thickness flux
641     C based on the average slope over the surface layer
642     C If we are on the edge of a "cliff" the surface layer at the
643     C centre of the grid point could be deeper than the U or V point.
644     C So, we ensure that we always take a sensible slope.
645     IF(kLow_U(i,j).LT.surfk(i,j)) THEN
646     kk=kLow_U(i,j)
647     hsurf = -rLowW(i,j,bi,bj)
648     ELSE
649     kk=surfk(i,j)
650     hsurf = -surfkz(i,j)
651     ENDIF
652     IF(kk.GT.0) THEN
653     IF(kk.EQ.Nr) THEN
654     tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
655     & *SlopeX(i,j,kk)/hsurf
656     ELSE
657     tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
658     & *( SlopeX(i,j,kk)-SlopeX(i,j,kk+1) )/hsurf
659     ENDIF
660     ELSE
661     tfluxX(i,j,k) = zeroRL
662     ENDIF
663    
664     IF(kLow_V(i,j).LT.surfk(i,j)) THEN
665     kk=kLow_V(i,j)
666     hsurf = -rLowS(i,j,bi,bj)
667     ELSE
668     kk=surfk(i,j)
669     hsurf = -surfkz(i,j)
670     ENDIF
671     IF(kk.GT.0) THEN
672     IF(kk.EQ.Nr) THEN
673     tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
674     & *SlopeY(i,j,kk)/hsurf
675     ELSE
676     tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
677     & *( SlopeY(i,j,kk)-SlopeY(i,j,kk+1) )/hsurf
678     ENDIF
679     ELSE
680     tfluxY(i,j,k) = zeroRL
681     ENDIF
682 m_bates 1.4 ENDIF
683 m_bates 1.1
684 m_bates 1.14 IF (GM_K3D_surfK) THEN
685     C Use a constant K in the surface layer.
686     KPV(i,j,k) = GM_K3D_constK
687 m_bates 1.4 ENDIF
688 m_bates 1.1 ENDIF
689 m_bates 1.14 ENDDO
690     ENDDO
691 m_bates 1.1 ENDDO
692 m_bates 1.14 ENDIF
693 m_bates 1.1
694     C Calculate gradq
695 m_bates 1.5 IF (GM_K3D_likeGM .OR. GM_K3D_beta_eq_0) THEN
696     C Ignore beta in the calculation of grad(q)
697 m_bates 1.2 DO k=1,Nr
698     DO j=1-Oly+1,sNy+Oly
699     DO i=1-Olx+1,sNx+Olx
700     gradqx(i,j,k) = maskW(i,j,k,bi,bj)*tfluxX(i,j,k)
701     gradqy(i,j,k) = maskS(i,j,k,bi,bj)*tfluxY(i,j,k)
702     ENDDO
703     ENDDO
704     ENDDO
705    
706     ELSE
707     C Do not ignore beta
708     DO k=1,Nr
709     DO j=1-Oly+1,sNy+Oly
710     DO i=1-Olx+1,sNx+Olx
711     gradqx(i,j,k) = maskW(i,j,k,bi,bj)*(dfdx(i,j)+tfluxX(i,j,k))
712     gradqy(i,j,k) = maskS(i,j,k,bi,bj)*(dfdy(i,j)+tfluxY(i,j,k))
713     ENDDO
714     ENDDO
715 m_bates 1.1 ENDDO
716 m_bates 1.2 ENDIF
717 m_bates 1.1
718     C ======================================
719     C Find Xi and the eddy induced velocities
720     C ======================================
721     C Find the buoyancy frequency at the west and south faces of a cell
722     C This is necessary to find the eigenvectors at those points
723     DO k=1,Nr
724     DO j=1-Oly+1,sNy+Oly
725     DO i=1-Olx+1,sNx+Olx
726     N2W(i,j,k) = maskW(i,j,k,bi,bj)
727     & *( N2(i,j,k)+N2(i-1,j,k) )
728     N2S(i,j,k) = maskS(i,j,k,bi,bj)
729     & *( N2(i,j,k)+N2(i,j-1,k) )
730     ENDDO
731     ENDDO
732     C This fudge is necessary to avoid division by zero in gmredi_calc_eigs.
733 jmc 1.7 C It does not affect the end result since it is in the overlap region.
734 m_bates 1.1 j=1-Oly
735     DO i=1-Olx,sNx+Olx
736     N2W(i,j,k) = GM_K3D_minN2
737     N2S(i,j,k) = GM_K3D_minN2
738     ENDDO
739     i=1-Olx
740     DO j=1-Oly,sNy+Oly
741     N2W(i,j,k) = GM_K3D_minN2
742     N2S(i,j,k) = GM_K3D_minN2
743     ENDDO
744     ENDDO
745    
746 m_bates 1.14 C If GM_K3D_likeGM=.TRUE., the diffusivity for the eddy transport is
747     C set to a constant equal to GM_K3D_constK.
748 m_bates 1.11 C If the diffusivity is constant the method here is the same as GM.
749 m_bates 1.14 C If GM_K3D_constRedi=.TRUE. K3D will be set equal to GM_K3D_constK later.
750 m_bates 1.2 IF(GM_K3D_likeGM) THEN
751 m_bates 1.3 DO k=1,Nr
752     DO j=1-Oly,sNy+Oly
753     DO i=1-Olx,sNx+Olx
754     KPV(i,j,k) = GM_K3D_constK
755     ENDDO
756 m_bates 1.2 ENDDO
757     ENDDO
758 m_bates 1.3 ENDIF
759 m_bates 1.2
760 m_bates 1.3 IF (.NOT. GM_K3D_smooth) THEN
761     C Do not expand K grad(q) => no smoothing
762     C May only be done with a constant K, otherwise the
763     C integral constraint is violated.
764 m_bates 1.2 DO k=1,Nr
765     DO j=1-Oly,sNy+Oly
766     DO i=1-Olx,sNx+Olx
767 m_bates 1.3 Xix(i,j,k) = -maskW(i,j,k,bi,bj)*KPV(i,j,k)*gradqx(i,j,k)
768     Xiy(i,j,k) = -maskS(i,j,k,bi,bj)*KPV(i,j,k)*gradqy(i,j,k)
769 m_bates 1.2 ENDDO
770     ENDDO
771     ENDDO
772    
773     ELSE
774 m_bates 1.3 C Expand K grad(q) in terms of baroclinic modes to smooth
775     C and satisfy the integral constraint
776 m_bates 1.2
777 m_bates 1.1 C Start with the X direction
778     C ------------------------------
779     C Calculate the eigenvectors at the West face of a cell
780 m_bates 1.4 IF (update_modes) THEN
781     CALL GMREDI_CALC_EIGS(
782     I iMin,iMax,jMin,jMax,bi,bj,N2W,myThid,
783     I kLow_U,maskW(:,:,:,bi,bj),
784     I hfacW(:,:,:,bi,bj),recip_hfacW(:,:,:,bi,bj),
785     I rLowW(:,:,bi,bj),GM_K3D_NModes,.FALSE.,
786     O dummy,modesW(:,:,:,:,bi,bj))
787     ENDIF
788 m_bates 1.1
789     C Calculate Xi_m at the west face of a cell
790     DO j=1-Oly,sNy+Oly
791     DO i=1-Olx,sNx+Olx
792     DO m=1,GM_K3D_NModes
793     XimX(m,i,j) = zeroRL
794     ENDDO
795     ENDDO
796     ENDDO
797     DO k=1,Nr
798     DO j=1-Oly,sNy+Oly
799     DO i=1-Olx,sNx+Olx
800     DO m=1,GM_K3D_NModes
801 m_bates 1.13 Kdqdx(i,j,k) = KPV(i,j,k)*gradqx(i,j,k)
802     XimX(m,i,j) = XimX(m,i,j)
803     & - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)
804     & *Kdqdx(i,j,k)*modesW(m,i,j,k,bi,bj)
805 m_bates 1.1 ENDDO
806     ENDDO
807     ENDDO
808     ENDDO
809    
810     C Calculate Xi in the X direction at the west face
811     DO k=1,Nr
812     DO j=1-Oly,sNy+Oly
813     DO i=1-Olx,sNx+Olx
814     Xix(i,j,k) = zeroRL
815     ENDDO
816     ENDDO
817     ENDDO
818     DO k=1,Nr
819     DO j=1-Oly,sNy+Oly
820     DO i=1-Olx,sNx+Olx
821     DO m=1,GM_K3D_NModes
822     Xix(i,j,k) = Xix(i,j,k)
823 m_bates 1.4 & + maskW(i,j,k,bi,bj)*XimX(m,i,j)*modesW(m,i,j,k,bi,bj)
824 m_bates 1.1 ENDDO
825     ENDDO
826     ENDDO
827     ENDDO
828    
829     C Now the Y direction
830     C ------------------------------
831     C Calculate the eigenvectors at the West face of a cell
832 m_bates 1.4 IF (update_modes) THEN
833     CALL GMREDI_CALC_EIGS(
834     I iMin,iMax,jMin,jMax,bi,bj,N2S,myThid,
835     I kLow_V,maskS(:,:,:,bi,bj),
836     I hfacS(:,:,:,bi,bj),recip_hfacS(:,:,:,bi,bj),
837     I rLowS(:,:,bi,bj), GM_K3D_NModes, .FALSE.,
838     O dummy,modesS(:,:,:,:,bi,bj))
839     ENDIF
840    
841 m_bates 1.1 DO j=1-Oly,sNy+Oly
842     DO i=1-Olx,sNx+Olx
843     DO m=1,GM_K3D_NModes
844     XimY(m,i,j) = zeroRL
845     ENDDO
846     ENDDO
847     ENDDO
848     DO k=1,Nr
849     DO j=1-Oly,sNy+Oly
850     DO i=1-Olx,sNx+Olx
851     DO m=1,GM_K3D_NModes
852 m_bates 1.13 Kdqdy(i,j,k) = KPV(i,j,k)*gradqy(i,j,k)
853 m_bates 1.3 XimY(m,i,j) = XimY(m,i,j)
854     & - drF(k)*hfacS(i,j,k,bi,bj)
855 m_bates 1.13 & *Kdqdy(i,j,k)*modesS(m,i,j,k,bi,bj)
856 m_bates 1.1 ENDDO
857     ENDDO
858     ENDDO
859     ENDDO
860    
861     C Calculate Xi for Y direction at the south face
862     DO k=1,Nr
863     DO j=1-Oly,sNy+Oly
864     DO i=1-Olx,sNx+Olx
865     Xiy(i,j,k) = zeroRL
866     ENDDO
867     ENDDO
868     ENDDO
869     DO k=1,Nr
870     DO j=1-Oly,sNy+Oly
871     DO i=1-Olx,sNx+Olx
872     DO m=1,GM_K3D_NModes
873     Xiy(i,j,k) = Xiy(i,j,k)
874 m_bates 1.4 & + maskS(i,j,k,bi,bj)*XimY(m,i,j)*modesS(m,i,j,k,bi,bj)
875 m_bates 1.1 ENDDO
876     ENDDO
877     ENDDO
878     ENDDO
879    
880 m_bates 1.11 C ENDIF (.NOT. GM_K3D_smooth)
881 m_bates 1.1 ENDIF
882    
883    
884     C Calculate the eddy induced velocity in the X direction at the west face
885     DO k=1,Nr
886     DO j=1-Oly+1,sNy+Oly
887     DO i=1-Olx+1,sNx+Olx
888     ustar(i,j,k) = -Xix(i,j,k)/coriU(i,j)
889     ENDDO
890     ENDDO
891     ENDDO
892    
893     C Calculate the eddy induced velocity in the Y direction at the south face
894     DO k=1,Nr
895     DO j=1-Oly+1,sNy+Oly
896     DO i=1-Olx+1,sNx+Olx
897     vstar(i,j,k) = -Xiy(i,j,k)/coriV(i,j)
898     ENDDO
899     ENDDO
900     ENDDO
901    
902     C ======================================
903     C Calculate the eddy induced overturning streamfunction
904     C ======================================
905     #ifdef GM_K3D_PASSIVE
906     k=Nr
907     DO j=1-Oly,sNy+Oly
908     DO i=1-Olx,sNx+Olx
909     psistar(i,j,Nr) = -hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
910     ENDDO
911     ENDDO
912     DO k=Nr-1,1,-1
913     DO j=1-Oly,sNy+Oly
914     DO i=1-Olx,sNx+Olx
915     psistar(i,j,k) = psistar(i,j,k+1)
916     & - hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
917     ENDDO
918     ENDDO
919     ENDDO
920    
921     #else
922    
923     IF (GM_AdvForm) THEN
924     k=Nr
925     DO j=1-Oly+1,sNy+1
926     DO i=1-Olx+1,sNx+1
927     GM_PsiX(i,j,k,bi,bj) = -hfacW(i,j,k,bi,bj)*drF(k)*ustar(i,j,k)
928     GM_PsiY(i,j,k,bi,bj) = -hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
929     ENDDO
930     ENDDO
931     DO k=Nr-1,1,-1
932     DO j=1-Oly+1,sNy+1
933     DO i=1-Olx+1,sNx+1
934     GM_PsiX(i,j,k,bi,bj) = GM_PsiX(i,j,k+1,bi,bj)
935     & - hfacW(i,j,k,bi,bj)*drF(k)*ustar(i,j,k)
936     GM_PsiY(i,j,k,bi,bj) = GM_PsiY(i,j,k+1,bi,bj)
937     & - hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
938     ENDDO
939     ENDDO
940     ENDDO
941    
942     ENDIF
943     #endif
944    
945     #ifdef ALLOW_DIAGNOSTICS
946     C Diagnostics
947     IF ( useDiagnostics ) THEN
948     CALL DIAGNOSTICS_FILL(K3D, 'GM_K3D ',0,Nr,0,1,1,myThid)
949 m_bates 1.14 CALL DIAGNOSTICS_FILL(KPV, 'GM_KPV ',0,Nr,0,1,1,myThid)
950 m_bates 1.1 CALL DIAGNOSTICS_FILL(urms, 'GM_URMS ',0,Nr,0,1,1,myThid)
951     CALL DIAGNOSTICS_FILL(Rdef, 'GM_RDEF ',0, 1,0,1,1,myThid)
952 m_bates 1.8 CALL DIAGNOSTICS_FILL(Rurms, 'GM_RURMS',0, 1,0,1,1,myThid)
953     CALL DIAGNOSTICS_FILL(RRhines,'GM_RRHNS',0, 1,0,1,1,myThid)
954 m_bates 1.1 CALL DIAGNOSTICS_FILL(Rmix, 'GM_RMIX ',0, 1,0,1,1,myThid)
955     CALL DIAGNOSTICS_FILL(supp, 'GM_SUPP ',0,Nr,0,1,1,myThid)
956     CALL DIAGNOSTICS_FILL(Xix, 'GM_Xix ',0,Nr,0,1,1,myThid)
957     CALL DIAGNOSTICS_FILL(Xiy, 'GM_Xiy ',0,Nr,0,1,1,myThid)
958     CALL DIAGNOSTICS_FILL(cDopp, 'GM_C ',0, 1,0,1,1,myThid)
959     CALL DIAGNOSTICS_FILL(Ubaro, 'GM_UBARO',0, 1,0,1,1,myThid)
960     CALL DIAGNOSTICS_FILL(eady, 'GM_EADY ',0, 1,0,1,1,myThid)
961     CALL DIAGNOSTICS_FILL(SlopeX, 'GM_Sx ',0,Nr,0,1,1,myThid)
962     CALL DIAGNOSTICS_FILL(SlopeY, 'GM_Sy ',0,Nr,0,1,1,myThid)
963     CALL DIAGNOSTICS_FILL(tfluxX, 'GM_TFLXX',0,Nr,0,1,1,myThid)
964     CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,0,1,1,myThid)
965     CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,0,1,1,myThid)
966     CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,0,1,1,myThid)
967 m_bates 1.13 CALL DIAGNOSTICS_FILL(Kdqdy, 'GM_Kdqdy',0,Nr,0,1,1,myThid)
968     CALL DIAGNOSTICS_FILL(Kdqdx, 'GM_Kdqdx',0,Nr,0,1,1,myThid)
969 m_bates 1.1 CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,0,1,1,myThid)
970     CALL DIAGNOSTICS_FILL(ustar, 'GM_USTAR',0,Nr,0,1,1,myThid)
971     CALL DIAGNOSTICS_FILL(vstar, 'GM_VSTAR',0,Nr,0,1,1,myThid)
972     CALL DIAGNOSTICS_FILL(umc, 'GM_UMC ',0,Nr,0,1,1,myThid)
973     CALL DIAGNOSTICS_FILL(ubar, 'GM_UBAR ',0,Nr,0,1,1,myThid)
974 m_bates 1.4 CALL DIAGNOSTICS_FILL(modesC(1,:,:,:,bi,bj),
975     & 'GM_MODEC',0,Nr,0,1,1,myThid)
976 m_bates 1.10 CALL DIAGNOSTICS_FILL(M4loc, 'GM_M4 ',0,Nr,0,1,1,myThid)
977     CALL DIAGNOSTICS_FILL(N2loc, 'GM_N2 ',0,Nr,0,1,1,myThid)
978 m_bates 1.13 CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,0,1,1,myThid)
979 m_bates 1.10 CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,0,1,1,myThid)
980    
981 m_bates 1.1 ENDIF
982     #endif
983    
984 m_bates 1.11 C For the Redi diffusivity, we set K3D to a constant if
985 m_bates 1.14 C GM_K3D_constRedi=.TRUE.
986     IF (GM_K3D_constRedi) THEN
987 m_bates 1.11 DO k=1,Nr
988     DO j=1-Oly,sNy+Oly
989     DO i=1-Olx,sNx+Olx
990     K3D(i,j,k,bi,bj) = GM_K3D_constK
991     ENDDO
992     ENDDO
993     ENDDO
994     ENDIF
995    
996 m_bates 1.14 #ifdef ALLOW_DIAGNOSTICS
997     IF ( useDiagnostics )
998     & CALL DIAGNOSTICS_FILL(K3D, 'GM_K3D_T',0,Nr,0,1,1,myThid)
999     #endif
1000    
1001 m_bates 1.1 #endif /* GM_K3D */
1002     RETURN
1003     END
  ViewVC Help
Powered by ViewVC 1.1.22