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