/[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.18 - (hide annotations) (download)
Tue May 20 11:42:28 2014 UTC (10 years ago) by m_bates
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint65
Changes since 1.17: +5 -5 lines 
for K3D parameterisation: fixed up minor bugs relating to finding the eigenvalues/vectors
1 m_bates 1.18 C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_k3d.F,v 1.17 2014/05/18 02:38:55 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 m_bates 1.16 C centreX,centreY :: used for calculating averages at centre of cell
154     C numerator,denominator :: of the renormalisation factor
155     C uInt :: column integral of u velocity (sum u*dz)
156     C vInt :: column integral of v velocity (sum v*dz)
157     C KdqdxInt :: column integral of K*dqdx (sum K*dqdx*dz)
158     C KdqdyInt :: column integral of K*dqdy (sum K*dqdy*dz)
159     C uKdqdyInt :: column integral of u*K*dqdy (sum u*K*dqdy*dz)
160     C vKdqdxInt :: column integral of v*K*dqdx (sum v*K*dqdx*dz)
161     C uXiyInt :: column integral of u*Xiy (sum u*Xiy*dz)
162     C vXixInt :: column integral of v*Xix (sum v*Xix*dz)
163     C Renorm :: renormalisation factor at the centre of a cell
164     C RenormU :: renormalisation factor at the western face of a cell
165     C RenormV :: renormalisation factor at the southern face of a cell
166 m_bates 1.15 _RL centreX, centreY
167     _RL numerator, denominator
168     _RL uInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
169     _RL vInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
170     _RL KdqdxInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
171     _RL KdqdyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
172     _RL uKdqdyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
173     _RL vKdqdxInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
174     _RL uXiyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
175     _RL vXixInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
176     _RL Renorm(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
177     _RL RenormU(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
178     _RL RenormV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
179    
180 m_bates 1.4 C gradqx :: Potential vorticity gradient in x direction
181     C gradqy :: Potential vorticity gradient in y direction
182     C XimX :: Vertical integral of phi_m*K*gradqx
183     C XimY :: Vertical integral of phi_m*K*gradqy
184     C cDopp :: Quasi-Doppler shifted long Rossby wave speed (m/s)
185     C umc :: ubar-c (m/s)
186     C eady :: Eady growth rate (1/s)
187     C urms :: the rms eddy velocity (m/s)
188     C supp :: The suppression factor
189     C ustar :: The eddy induced velocity in the x direction
190     C vstar :: The eddy induced velocity in the y direction
191     C Xix :: Xi in the x direction (m/s**2)
192     C Xiy :: Xi in the y direction (m/s**2)
193 m_bates 1.1 _RL gradqx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
194     _RL gradqy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
195     _RL XimX(GM_K3D_NModes,1-Olx:sNx+Olx,1-Oly:sNy+Oly)
196     _RL XimY(GM_K3D_NModes,1-Olx:sNx+Olx,1-Oly:sNy+Oly)
197 m_bates 1.4 _RL cDopp(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
198 m_bates 1.1 _RL umc( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
199     _RL eady( 1-Olx:sNx+Olx,1-Oly:sNy+Oly)
200     _RL urms( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
201     _RL supp( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
202     _RL ustar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
203     _RL vstar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
204     _RL Xix( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
205     _RL Xiy( 1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
206 m_bates 1.4 #ifdef GM_K3D_PASSIVE
207     C psistar :: eddy induced streamfunction in the y direction
208 m_bates 1.1 _RL psistar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,1:Nr)
209 m_bates 1.4 #endif
210 m_bates 1.1
211    
212     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
213    
214     C ======================================
215     C Initialise some variables
216     C ======================================
217     small = TINY(zeroRL)
218 m_bates 1.4 update_modes=.FALSE.
219     IF ( DIFFERENT_MULTIPLE(GM_K3D_vecFreq,myTime,deltaTClock) )
220     & update_modes=.TRUE.
221    
222 m_bates 1.1 DO j=1-Oly,sNy+Oly
223     DO i=1-Olx,sNx+Olx
224     kLow_C(i,j) = kLowC(i,j,bi,bj)
225     ENDDO
226     ENDDO
227     DO j=1-Oly,sNy+Oly
228     DO i=1-Olx+1,sNx+Olx
229     kLow_U(i,j) = MIN( kLow_C(i,j), kLow_C(i-1,j) )
230     ENDDO
231     ENDDO
232     DO j=1-Oly+1,sNy+Oly
233     DO i=1-Olx,sNx+Olx
234     kLow_V(i,j) = MIN( kLow_C(i,j), kLow_C(i,j-1) )
235     ENDDO
236     ENDDO
237    
238 m_bates 1.18 C Dummy values for the edges. This does not affect the results
239     C but avoids problems when solving for the eigenvalues.
240 m_bates 1.1 i=1-Olx
241     DO j=1-Oly,sNy+Oly
242 m_bates 1.18 kLow_U(i,j) = 0
243 m_bates 1.1 ENDDO
244     j=1-Oly
245     DO i=1-Olx,sNx+Olx
246 m_bates 1.18 kLow_V(i,j) = 0
247 m_bates 1.1 ENDDO
248    
249     g_reciprho_sq = (gravity*recip_rhoConst)**2
250     C Gradient of Coriolis
251     DO j=1-Oly+1,sNy+Oly
252     DO i=1-Olx+1,sNx+Olx
253     dfdx(i,j) = ( fCori(i,j,bi,bj)-fCori(i-1,j,bi,bj) )
254     & *recip_dxC(i,j,bi,bj)
255     dfdy(i,j) = ( fCori(i,j,bi,bj)-fCori(i,j-1,bi,bj) )
256     & *recip_dyC(i,j,bi,bj)
257     ENDDO
258     ENDDO
259    
260     C Coriolis at C points enforcing a minimum value so
261     C that it is defined at the equator
262     DO j=1-Oly,sNy+Oly
263     DO i=1-Olx,sNx+Olx
264     cori(i,j) = SIGN( MAX( ABS(fCori(i,j,bi,bj)),GM_K3D_minCori ),
265     & fCori(i,j,bi,bj) )
266     ENDDO
267     ENDDO
268     C Coriolis at U and V points
269 m_bates 1.17 DO j=1-Oly,sNy+Oly
270 m_bates 1.1 DO i=1-Olx+1,sNx+Olx
271     C Limited so that the inverse is defined at the equator
272     coriU(i,j) = op5*( cori(i,j)+cori(i-1,j) )
273     coriU(i,j) = SIGN( MAX( ABS(coriU(i,j)),GM_K3D_minCori ),
274     & coriU(i,j) )
275    
276 m_bates 1.17 C Not limited
277     fCoriU(i,j) = op5*( fCori(i,j,bi,bj)+fCori(i-1,j,bi,bj) )
278     ENDDO
279     ENDDO
280     DO j=1-Oly+1,sNy+Oly
281     DO i=1-Olx,sNx+Olx
282     C Limited so that the inverse is defined at the equator
283 m_bates 1.1 coriV(i,j) = op5*( cori(i,j)+cori(i,j-1) )
284     coriV(i,j) = SIGN( MAX( ABS(coriV(i,j)),GM_K3D_minCori ),
285     & coriV(i,j) )
286    
287     C Not limited
288     fCoriV(i,j) = op5*( fCori(i,j,bi,bj)+fCori(i,j-1,bi,bj) )
289     ENDDO
290     ENDDO
291     DO j=1-Oly,sNy+Oly
292     DO i=1-Olx,sNx+Olx
293 m_bates 1.6 gradf(i,j) = recip_rSphere*fCoriCos(i,j,bi,bj)
294 m_bates 1.1 ENDDO
295     ENDDO
296 m_bates 1.17 C Some dummy values at the edges
297     i=1-Olx
298     DO j=1-Oly,sNy+Oly
299     coriU(i,j)=cori(i,j)
300     fCoriU(i,j)=fCori(i,j,bi,bj)
301     ENDDO
302     j=1-Oly
303     DO i=1-Olx,sNx+Olx
304     coriV(i,j)=cori(i,j)
305     fCoriV(i,j)=fCori(i,j,bi,bj)
306     ENDDO
307 m_bates 1.1
308     C Zeroing some cumulative fields
309     DO j=1-Oly,sNy+Oly
310     DO i=1-Olx,sNx+Olx
311 m_bates 1.10 eady(i,j) = zeroRL
312     BVint(i,j) = zeroRL
313     Ubaro(i,j) = zeroRL
314     deltaH(i,j) = zeroRL
315     ENDDO
316     ENDDO
317     DO k=1,Nr
318     DO j=1-Oly,sNy+Oly
319     DO i=1-Olx,sNx+Olx
320     slopeC(i,j,k)=zeroRL
321     ENDDO
322 m_bates 1.1 ENDDO
323     ENDDO
324    
325 m_bates 1.17 C initialise remaining 2d variables
326     DO j=1-Oly,sNy+Oly
327     DO i=1-Olx,sNx+Olx
328     dfdy(i,j)=zeroRL
329     dfdy(i,j)=zeroRL
330     Rurms(i,j)=zeroRL
331     RRhines(i,j)=zeroRL
332     Rmix(i,j)=zeroRL
333     ENDDO
334     ENDDO
335     C initialise remaining 3d variables
336     DO k=1,Nr
337     DO j=1-Oly,sNy+Oly
338     DO i=1-Olx,sNx+Olx
339     N2loc(i,j,k)=GM_K3D_minN2
340     N2W(i,j,k) = GM_K3D_minN2
341     N2S(i,j,k) = GM_K3D_minN2
342     M4loc(i,j,k)=zeroRL
343     M4onN2(i,j,k)=zeroRL
344     urms(i,j,k)=zeroRL
345     SlopeX(i,j,k)=zeroRL
346     SlopeY(i,j,k)=zeroRL
347     dSigmaDx(i,j,k)=zeroRL
348     dSigmaDy(i,j,k)=zeroRL
349     gradqx(i,j,k)=zeroRL
350     gradqy(i,j,k)=zeroRL
351     ENDDO
352     ENDDO
353     ENDDO
354    
355 m_bates 1.1 C Find the zonal velocity at the cell centre
356     C The logicals here are, in order: 1/ go from grid to north/east directions
357     C 2/ go from C to A grid and 3/ apply the mask
358 m_bates 1.9 #ifdef ALLOW_EDDYPSI
359     IF (GM_InMomAsStress) THEN
360     CALL rotate_uv2en_rl(uMean, vMean, ubar, vbar, .TRUE., .TRUE.,
361     & .TRUE.,Nr,mythid)
362     ELSE
363     #endif
364     CALL rotate_uv2en_rl(uVel, vVel, ubar, vbar, .TRUE., .TRUE.,
365     & .TRUE.,Nr,mythid)
366     #ifdef ALLOW_EDDYPSI
367     ENDIF
368     #endif
369 m_bates 1.1
370     C Square of the buoyancy frequency at the top of a grid cell
371     DO k=2,Nr
372     DO j=1-Oly,sNy+Oly
373     DO i=1-Olx,sNx+Olx
374     N2(i,j,k) = -gravity*recip_rhoConst*sigmaR(i,j,k)
375     ENDDO
376     ENDDO
377     ENDDO
378     C N2(k=1) is always zero
379     k=1
380     DO j=1-Oly,sNy+Oly
381     DO i=1-Olx,sNx+Olx
382 m_bates 1.17 N2(i,j,k) = zeroRL
383     N(i,j,k) = zeroRL
384 m_bates 1.1 ENDDO
385     ENDDO
386     C Enforce a minimum N2
387     DO k=2,Nr
388     DO j=1-Oly,sNy+Oly
389     DO i=1-Olx,sNx+Olx
390     IF (N2(i,j,k).LT.GM_K3D_minN2) N2(i,j,k)=GM_K3D_minN2
391     N(i,j,k) = SQRT(N2(i,j,k))
392     ENDDO
393     ENDDO
394     ENDDO
395     C Calculate the minimum drho/dz
396     maxDRhoDz = -rhoConst*GM_K3D_minN2/gravity
397    
398     C Calculate the barotropic velocity by vertically integrating
399     C and the dividing by the depth of the water column
400     C Note that Ubaro is on the U grid.
401     DO k=1,Nr
402     DO j=1-Oly,sNy+Oly
403     DO i=1-Olx,sNx+Olx
404     Ubaro(i,j) = Ubaro(i,j) +
405     & maskW(i,j,k,bi,bj)*drF(k)*hfacC(i,j,k,bi,bj)
406     & *ubar(i,j,k,bi,bj)
407     ENDDO
408     ENDDO
409     ENDDO
410     DO j=1-Oly,sNy+Oly
411     DO i=1-Olx,sNx+Olx
412     IF (kLow_C(i,j).GT.0) THEN
413     C The minus sign is because r_Low<0
414     Ubaro(i,j) = -Ubaro(i,j)/r_Low(i,j,bi,bj)
415     ENDIF
416     ENDDO
417     ENDDO
418    
419     C Integrate the buoyancy frequency vertically using the trapezoidal method.
420     DO k=1,Nr
421     DO j=1-Oly,sNy+Oly
422     DO i=1-Olx,sNx+Olx
423     IF (k.LT.kLow_C(i,j)) THEN
424     BVint(i,j) = BVint(i,j) + hFacC(i,j,k,bi,bj)*drF(k)
425     & *(N(i,j,k)+N(i,j,k+1))
426     ELSEIF (k.EQ.kLow_C(i,j)) THEN
427     C Assume that N(z=-H)=0
428     BVint(i,j) = BVint(i,j) + hFacC(i,j,k,bi,bj)*drF(k)*N(i,j,k)
429     ENDIF
430     ENDDO
431     ENDDO
432     ENDDO
433     DO j=1-Oly,sNy+Oly
434     DO i=1-Olx,sNx+Olx
435     BVint(i,j) = op5*BVint(i,j)
436     ENDDO
437     ENDDO
438    
439     C Calculate the eigenvalues and eigenvectors
440 m_bates 1.4 IF (update_modes) THEN
441     CALL GMREDI_CALC_EIGS(
442     I iMin,iMax,jMin,jMax,bi,bj,N2,myThid,
443     I kLow_C, maskC(:,:,:,bi,bj),
444     I hfacC(:,:,:,bi,bj), recip_hfacC(:,:,:,bi,bj),
445     I R_Low(:,:,bi,bj), 1, .TRUE.,
446     O Rmid, modesC(:,:,:,:,bi,bj))
447    
448     C Calculate the Rossby Radius
449     DO j=1-Oly+1,sNy+Oly
450     DO i=1-Olx+1,sNx+Olx
451     Req = SQRT(BVint(i,j)/(2*pi*gradf(i,j)))
452     Rdef(i,j,bi,bj) = MIN(Rmid(i,j),Req)
453     ENDDO
454     ENDDO
455     ENDIF
456 m_bates 1.1
457     C Average dsigma/dx and dsigma/dy onto the centre points
458    
459     DO k=1,Nr
460     DO j=1-Oly,sNy+Oly-1
461     DO i=1-Olx,sNx+Olx-1
462     dSigmaDx(i,j,k) = op5*(sigmaX(i,j,k)+sigmaX(i+1,j,k))
463     dSigmaDy(i,j,k) = op5*(sigmaY(i,j,k)+sigmaY(i,j+1,k))
464     ENDDO
465     ENDDO
466     ENDDO
467    
468     C ===============================
469     C Calculate the Eady growth rate
470     C ===============================
471     DO k=1,Nr
472    
473 m_bates 1.10 DO j=1-Oly,sNy+Oly-1
474     DO i=1-Olx,sNx+Olx-1
475     M4loc(i,j,k) = g_reciprho_sq*( dSigmaDx(i,j,k)**2
476     & +dSigmaDy(i,j,k)**2 )
477     IF (k.NE.kLow_C(i,j)) THEN
478     N2loc(i,j,k) = op5*(N2(i,j,k)+N2(i,j,k+1))
479     ELSE
480     N2loc(i,j,k) = op5*N2(i,j,k)
481     ENDIF
482     ENDDO
483     ENDDO
484 m_bates 1.1 C The bottom of the grid cell is shallower than the top
485     C integration level, so, advance the depth.
486 m_bates 1.10 IF (-rF(k+1) .LE. GM_K3D_EadyMinDepth) CYCLE
487 m_bates 1.1
488 jmc 1.7 C Do not bother going any deeper since the top of the
489 m_bates 1.1 C cell is deeper than the bottom integration level
490     IF (-rF(k).GE.GM_K3D_EadyMaxDepth) EXIT
491    
492     C We are in the integration depth range
493     DO j=1-Oly,sNy+Oly-1
494     DO i=1-Olx,sNx+Olx-1
495 m_bates 1.10 IF ( (kLow_C(i,j).GE.k) .AND.
496     & (-hMixLayer(i,j,bi,bj).LE.-rC(k)) ) THEN
497 m_bates 1.1
498 m_bates 1.12 slopeC(i,j,k) = SQRT(M4loc(i,j,k))/N2loc(i,j,k)
499 m_bates 1.1 C Limit the slope. Note, this is not all the Eady calculations.
500 m_bates 1.13 IF (slopeC(i,j,k).LE.GM_maxSlope) THEN
501     M4onN2(i,j,k) = M4loc(i,j,k)/N2loc(i,j,k)
502 m_bates 1.1 ELSE
503 m_bates 1.13 slopeC(i,j,k) = GM_maxslope
504     M4onN2(i,j,k) = SQRT(M4loc(i,j,k))*GM_maxslope
505 m_bates 1.1 ENDIF
506 m_bates 1.13 eady(i,j) = eady(i,j)
507     & + hfacC(i,j,k,bi,bj)*drF(k)*M4onN2(i,j,k)
508 m_bates 1.10 deltaH(i,j) = deltaH(i,j) + drF(k)
509 m_bates 1.1 ENDIF
510     ENDDO
511     ENDDO
512     ENDDO
513    
514     DO j=1-Oly,sNy+Oly
515     DO i=1-Olx,sNx+Olx
516 m_bates 1.10 C If the minimum depth for the integration is deeper than the ocean
517     C bottom OR the mixed layer is deeper than the maximum depth of
518     C integration, we set the Eady growth rate to something small
519     C to avoid floating point exceptions.
520     C Later, these areas will be given a small diffusivity.
521     IF (deltaH(i,j).EQ.zeroRL) THEN
522 m_bates 1.1 eady(i,j) = small
523    
524 m_bates 1.10 C Otherwise, divide over the integration and take the square root
525     C to actually find the Eady growth rate.
526 m_bates 1.1 ELSE
527 m_bates 1.10 eady(i,j) = SQRT(eady(i,j)/deltaH(i,j))
528 m_bates 1.1
529     ENDIF
530    
531     ENDDO
532     ENDDO
533    
534     C ======================================
535     C Calculate the diffusivity
536     C ======================================
537     DO j=1-Oly+1,sNy+Oly
538     DO i=1-Olx+1,sNx+Olx-1
539     C Calculate the Visbeck velocity
540 m_bates 1.13 Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_K3D_Rmax)
541 m_bates 1.8 urms(i,j,1) = GM_K3D_Lambda*eady(i,j)*Rurms(i,j)
542 m_bates 1.1 C Set the bottom urms to zero
543     k=kLow_C(i,j)
544     IF (k.GT.0) urms(i,j,k) = 0.0
545    
546     C Calculate the Rhines scale
547     RRhines(i,j) = SQRT(urms(i,j,1)/gradf(i,j))
548    
549     C Calculate the estimated length scale
550 m_bates 1.4 Rmix(i,j) = MIN(Rdef(i,j,bi,bj), RRhines(i,j))
551 m_bates 1.13 Rmix(i,j) = MAX(Rmix(i,j),GM_K3D_Rmin)
552 m_bates 1.1
553     C Calculate the Doppler shifted long Rossby wave speed
554     C Ubaro is on the U grid so we must average onto the M grid.
555     cDopp(i,j) = op5*( Ubaro(i,j)+Ubaro(i+1,j) )
556 m_bates 1.4 & - gradf(i,j)*Rdef(i,j,bi,bj)*Rdef(i,j,bi,bj)
557 m_bates 1.1 C Limit the wave speed to the namelist variable GM_K3D_maxC
558     IF (ABS(cDopp(i,j)).GT.GM_K3D_maxC) THEN
559     cDopp(i,j) = MAX(GM_K3D_maxC, cDopp(i,j))
560     ENDIF
561    
562     ENDDO
563     ENDDO
564    
565     C Project the surface urms to the subsurface using the first baroclinic mode
566 m_bates 1.4 CALL GMREDI_CALC_URMS(
567     I iMin,iMax,jMin,jMax,bi,bj,N2,myThid,
568     U urms)
569 m_bates 1.1
570     C Calculate the diffusivity (on the mass grid)
571     DO k=1,Nr
572     DO j=1-Oly,sNy+Oly
573     DO i=1-Olx,sNx+Olx
574     IF (k.LE.kLow_C(i,j)) THEN
575 m_bates 1.10 IF (deltaH(i,j).EQ.zeroRL) THEN
576 m_bates 1.1 K3D(i,j,k,bi,bj) = GM_K3D_smallK
577     ELSE
578     IF (urms(i,j,k).EQ.0.0) THEN
579     K3D(i,j,k,bi,bj) = GM_K3D_smallK
580     ELSE
581 m_bates 1.13 umc(i,j,k) =ubar(i,j,k,bi,bj) - cDopp(i,j)
582     supp(i,j,k)=1/(1+GM_K3D_b1*umc(i,j,k)**2/urms(i,j,1)**2)
583     C 2*Rmix gives the diameter
584 m_bates 1.1 K3D(i,j,k,bi,bj) = GM_K3D_gamma*urms(i,j,k)
585 m_bates 1.13 & *2*Rmix(i,j)*supp(i,j,k)
586 m_bates 1.1 ENDIF
587    
588     C Enforce lower and upper bounds on the diffusivity
589     IF (K3D(i,j,k,bi,bj).LT.GM_K3D_smallK)
590     & K3D(i,j,k,bi,bj) = GM_K3D_smallK
591     IF (K3D(i,j,k,bi,bj).GT.GM_maxK3D)
592     & K3D(i,j,k,bi,bj) = GM_maxK3D
593     ENDIF
594     ENDIF
595     ENDDO
596     ENDDO
597     ENDDO
598    
599     C ======================================
600     C Find the PV gradient
601     C ======================================
602 m_bates 1.3 C Calculate the surface layer thickness.
603     C Use hMixLayer (calculated in model/src/calc_oce_mxlayer)
604     C for the mixed layer depth.
605 m_bates 1.1
606 m_bates 1.3 C Enforce a minimum surface layer depth
607 m_bates 1.1 DO j=1-Oly,sNy+Oly
608     DO i=1-Olx,sNx+Olx
609 m_bates 1.3 surfkz(i,j) = MIN(-GM_K3D_surfMinDepth,-hMixLayer(i,j,bi,bj))
610     surfkz(i,j) = MAX(surfkz(i,j),R_low(i,j,bi,bj))
611     IF(maskC(i,j,1,bi,bj).EQ.0.0) surfkz(i,j)=0.0
612     surfk(i,j) = 0
613 m_bates 1.1 ENDDO
614     ENDDO
615 m_bates 1.4 DO k=1,Nr
616 m_bates 1.1 DO j=1-Oly,sNy+Oly
617     DO i=1-Olx,sNx+Olx
618 m_bates 1.3 IF (rF(k).GT.surfkz(i,j) .AND. surfkz(i,j).GE.rF(k+1))
619     & surfk(i,j) = k
620 m_bates 1.1 ENDDO
621     ENDDO
622     ENDDO
623 m_bates 1.3
624 m_bates 1.1 C Calculate the isopycnal slope
625     DO j=1-Oly,sNy+Oly-1
626     DO i=1-Olx,sNx+Olx-1
627     SlopeX(i,j,1) = zeroRL
628     SlopeY(i,j,1) = zeroRL
629     ENDDO
630     ENDDO
631     DO k=2,Nr
632     DO j=1-Oly+1,sNy+Oly
633     DO i=1-Olx+1,sNx+Olx
634     IF(surfk(i,j).GE.kLowC(i,j,bi,bj)) THEN
635     C If the surface layer is thinner than the water column
636     C the set the slope to zero to avoid problems.
637     SlopeX(i,j,k) = zeroRL
638     SlopeY(i,j,k) = zeroRL
639    
640     ELSE
641     C Calculate the zonal slope at the western cell face (U grid)
642 m_bates 1.4 sigz = MIN( op5*(sigmaR(i,j,k)+sigmaR(i-1,j,k)), maxDRhoDz )
643 m_bates 1.1 sigx = op5*( sigmaX(i,j,k)+sigmaX(i,j,k-1) )
644     slope = sigx/sigz
645     IF(ABS(slope).GT.GM_maxSlope)
646     & slope = SIGN(GM_maxSlope,slope)
647     SlopeX(i,j,k)=-maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)*slope
648    
649     C Calculate the meridional slope at the southern cell face (V grid)
650 m_bates 1.4 sigz = MIN( op5*(sigmaR(i,j,k)+sigmaR(i,j-1,k)), maxDRhoDz )
651 m_bates 1.1 sigy = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )
652     slope = sigy/sigz
653     IF(ABS(slope).GT.GM_maxSlope)
654     & slope = SIGN(GM_maxSlope,slope)
655     SlopeY(i,j,k)=-maskS(i,j,k-1,bi,bj)*maskS(i,j,k,bi,bj)*slope
656     ENDIF
657     ENDDO
658     ENDDO
659     ENDDO
660    
661 m_bates 1.14 C Calculate the thickness flux and diffusivity. These may be altered later
662     C depending on namelist options.
663 m_bates 1.1 C Enforce a zero slope bottom boundary condition for the bottom most cells (k=Nr)
664     k=Nr
665     DO j=1-Oly,sNy+Oly
666     DO i=1-Olx,sNx+Olx
667     C Zonal thickness flux at the western cell face
668     tfluxX(i,j,k) = -fCoriU(i,j)*SlopeX(i,j,k)
669     & *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
670     C Meridional thickness flux at the southern cell face
671     tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)
672     & *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
673 m_bates 1.14
674     C Use the interior diffusivity. Note: if we are using a
675     C constant diffusivity KPV is overwritten later
676     KPV(i,j,k) = K3D(i,j,k,bi,bj)
677    
678 m_bates 1.1 ENDDO
679     ENDDO
680    
681 m_bates 1.14 C Calculate the thickness flux and diffusivity and for other cells (k<Nr)
682 m_bates 1.1 DO k=Nr-1,1,-1
683 m_bates 1.14 DO j=1-Oly,sNy+Oly
684     DO i=1-Olx,sNx+Olx
685     C Zonal thickness flux at the western cell face
686     tfluxX(i,j,k)=-fCoriU(i,j)*(SlopeX(i,j,k)-SlopeX(i,j,k+1))
687     & *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
688     & *maskW(i,j,k,bi,bj)
689    
690     C Meridional thickness flux at the southern cell face
691     tfluxY(i,j,k)=-fCoriV(i,j)*(SlopeY(i,j,k)-SlopeY(i,j,k+1))
692     & *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
693     & *maskS(i,j,k,bi,bj)
694    
695     C Use the interior diffusivity. Note: if we are using a
696     C constant diffusivity KPV is overwritten later
697     KPV(i,j,k) = K3D(i,j,k,bi,bj)
698     ENDDO
699     ENDDO
700     ENDDO
701    
702    
703     C Special treatment for the surface layer (if necessary), which overwrites
704     C values in the previous loops.
705     IF (GM_K3D_ThickSheet .OR. GM_K3D_surfK) THEN
706     DO k=Nr-1,1,-1
707     DO j=1-Oly,sNy+Oly
708     DO i=1-Olx,sNx+Olx
709     IF(k.LE.surfk(i,j)) THEN
710     C We are in the surface layer. Change the thickness flux
711     C and diffusivity as necessary.
712    
713     IF (GM_K3D_ThickSheet) THEN
714     C We are in the surface layer, so set the thickness flux
715     C based on the average slope over the surface layer
716     C If we are on the edge of a "cliff" the surface layer at the
717     C centre of the grid point could be deeper than the U or V point.
718     C So, we ensure that we always take a sensible slope.
719     IF(kLow_U(i,j).LT.surfk(i,j)) THEN
720     kk=kLow_U(i,j)
721     hsurf = -rLowW(i,j,bi,bj)
722     ELSE
723     kk=surfk(i,j)
724     hsurf = -surfkz(i,j)
725     ENDIF
726     IF(kk.GT.0) THEN
727     IF(kk.EQ.Nr) THEN
728     tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
729     & *SlopeX(i,j,kk)/hsurf
730     ELSE
731     tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
732     & *( SlopeX(i,j,kk)-SlopeX(i,j,kk+1) )/hsurf
733     ENDIF
734     ELSE
735     tfluxX(i,j,k) = zeroRL
736     ENDIF
737    
738     IF(kLow_V(i,j).LT.surfk(i,j)) THEN
739     kk=kLow_V(i,j)
740     hsurf = -rLowS(i,j,bi,bj)
741     ELSE
742     kk=surfk(i,j)
743     hsurf = -surfkz(i,j)
744     ENDIF
745     IF(kk.GT.0) THEN
746     IF(kk.EQ.Nr) THEN
747     tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
748     & *SlopeY(i,j,kk)/hsurf
749     ELSE
750     tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
751     & *( SlopeY(i,j,kk)-SlopeY(i,j,kk+1) )/hsurf
752     ENDIF
753     ELSE
754     tfluxY(i,j,k) = zeroRL
755     ENDIF
756 m_bates 1.4 ENDIF
757 m_bates 1.1
758 m_bates 1.14 IF (GM_K3D_surfK) THEN
759     C Use a constant K in the surface layer.
760     KPV(i,j,k) = GM_K3D_constK
761 m_bates 1.4 ENDIF
762 m_bates 1.1 ENDIF
763 m_bates 1.14 ENDDO
764     ENDDO
765 m_bates 1.1 ENDDO
766 m_bates 1.14 ENDIF
767 m_bates 1.1
768     C Calculate gradq
769 m_bates 1.5 IF (GM_K3D_likeGM .OR. GM_K3D_beta_eq_0) THEN
770     C Ignore beta in the calculation of grad(q)
771 m_bates 1.2 DO k=1,Nr
772     DO j=1-Oly+1,sNy+Oly
773     DO i=1-Olx+1,sNx+Olx
774     gradqx(i,j,k) = maskW(i,j,k,bi,bj)*tfluxX(i,j,k)
775     gradqy(i,j,k) = maskS(i,j,k,bi,bj)*tfluxY(i,j,k)
776     ENDDO
777     ENDDO
778     ENDDO
779    
780     ELSE
781     C Do not ignore beta
782     DO k=1,Nr
783     DO j=1-Oly+1,sNy+Oly
784     DO i=1-Olx+1,sNx+Olx
785     gradqx(i,j,k) = maskW(i,j,k,bi,bj)*(dfdx(i,j)+tfluxX(i,j,k))
786     gradqy(i,j,k) = maskS(i,j,k,bi,bj)*(dfdy(i,j)+tfluxY(i,j,k))
787     ENDDO
788     ENDDO
789 m_bates 1.1 ENDDO
790 m_bates 1.2 ENDIF
791 m_bates 1.1
792     C ======================================
793     C Find Xi and the eddy induced velocities
794     C ======================================
795     C Find the buoyancy frequency at the west and south faces of a cell
796     C This is necessary to find the eigenvectors at those points
797     DO k=1,Nr
798     DO j=1-Oly+1,sNy+Oly
799     DO i=1-Olx+1,sNx+Olx
800     N2W(i,j,k) = maskW(i,j,k,bi,bj)
801     & *( N2(i,j,k)+N2(i-1,j,k) )
802     N2S(i,j,k) = maskS(i,j,k,bi,bj)
803     & *( N2(i,j,k)+N2(i,j-1,k) )
804     ENDDO
805     ENDDO
806     ENDDO
807    
808 m_bates 1.14 C If GM_K3D_likeGM=.TRUE., the diffusivity for the eddy transport is
809     C set to a constant equal to GM_K3D_constK.
810 m_bates 1.11 C If the diffusivity is constant the method here is the same as GM.
811 m_bates 1.14 C If GM_K3D_constRedi=.TRUE. K3D will be set equal to GM_K3D_constK later.
812 m_bates 1.2 IF(GM_K3D_likeGM) THEN
813 m_bates 1.3 DO k=1,Nr
814     DO j=1-Oly,sNy+Oly
815     DO i=1-Olx,sNx+Olx
816     KPV(i,j,k) = GM_K3D_constK
817     ENDDO
818 m_bates 1.2 ENDDO
819     ENDDO
820 m_bates 1.3 ENDIF
821 m_bates 1.2
822 m_bates 1.3 IF (.NOT. GM_K3D_smooth) THEN
823     C Do not expand K grad(q) => no smoothing
824     C May only be done with a constant K, otherwise the
825     C integral constraint is violated.
826 m_bates 1.2 DO k=1,Nr
827     DO j=1-Oly,sNy+Oly
828     DO i=1-Olx,sNx+Olx
829 m_bates 1.3 Xix(i,j,k) = -maskW(i,j,k,bi,bj)*KPV(i,j,k)*gradqx(i,j,k)
830     Xiy(i,j,k) = -maskS(i,j,k,bi,bj)*KPV(i,j,k)*gradqy(i,j,k)
831 m_bates 1.2 ENDDO
832     ENDDO
833     ENDDO
834    
835     ELSE
836 m_bates 1.3 C Expand K grad(q) in terms of baroclinic modes to smooth
837     C and satisfy the integral constraint
838 m_bates 1.2
839 m_bates 1.1 C Start with the X direction
840     C ------------------------------
841     C Calculate the eigenvectors at the West face of a cell
842 m_bates 1.4 IF (update_modes) THEN
843     CALL GMREDI_CALC_EIGS(
844     I iMin,iMax,jMin,jMax,bi,bj,N2W,myThid,
845     I kLow_U,maskW(:,:,:,bi,bj),
846     I hfacW(:,:,:,bi,bj),recip_hfacW(:,:,:,bi,bj),
847     I rLowW(:,:,bi,bj),GM_K3D_NModes,.FALSE.,
848     O dummy,modesW(:,:,:,:,bi,bj))
849     ENDIF
850 m_bates 1.1
851     C Calculate Xi_m at the west face of a cell
852     DO j=1-Oly,sNy+Oly
853     DO i=1-Olx,sNx+Olx
854     DO m=1,GM_K3D_NModes
855     XimX(m,i,j) = zeroRL
856     ENDDO
857     ENDDO
858     ENDDO
859     DO k=1,Nr
860     DO j=1-Oly,sNy+Oly
861     DO i=1-Olx,sNx+Olx
862     DO m=1,GM_K3D_NModes
863 m_bates 1.13 Kdqdx(i,j,k) = KPV(i,j,k)*gradqx(i,j,k)
864     XimX(m,i,j) = XimX(m,i,j)
865     & - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)
866     & *Kdqdx(i,j,k)*modesW(m,i,j,k,bi,bj)
867 m_bates 1.1 ENDDO
868     ENDDO
869     ENDDO
870     ENDDO
871    
872     C Calculate Xi in the X direction at the west face
873     DO k=1,Nr
874     DO j=1-Oly,sNy+Oly
875     DO i=1-Olx,sNx+Olx
876     Xix(i,j,k) = zeroRL
877     ENDDO
878     ENDDO
879     ENDDO
880     DO k=1,Nr
881     DO j=1-Oly,sNy+Oly
882     DO i=1-Olx,sNx+Olx
883     DO m=1,GM_K3D_NModes
884     Xix(i,j,k) = Xix(i,j,k)
885 m_bates 1.4 & + maskW(i,j,k,bi,bj)*XimX(m,i,j)*modesW(m,i,j,k,bi,bj)
886 m_bates 1.1 ENDDO
887     ENDDO
888     ENDDO
889     ENDDO
890    
891     C Now the Y direction
892     C ------------------------------
893     C Calculate the eigenvectors at the West face of a cell
894 m_bates 1.4 IF (update_modes) THEN
895     CALL GMREDI_CALC_EIGS(
896     I iMin,iMax,jMin,jMax,bi,bj,N2S,myThid,
897     I kLow_V,maskS(:,:,:,bi,bj),
898     I hfacS(:,:,:,bi,bj),recip_hfacS(:,:,:,bi,bj),
899     I rLowS(:,:,bi,bj), GM_K3D_NModes, .FALSE.,
900     O dummy,modesS(:,:,:,:,bi,bj))
901     ENDIF
902    
903 m_bates 1.1 DO j=1-Oly,sNy+Oly
904     DO i=1-Olx,sNx+Olx
905     DO m=1,GM_K3D_NModes
906     XimY(m,i,j) = zeroRL
907     ENDDO
908     ENDDO
909     ENDDO
910     DO k=1,Nr
911     DO j=1-Oly,sNy+Oly
912     DO i=1-Olx,sNx+Olx
913     DO m=1,GM_K3D_NModes
914 m_bates 1.13 Kdqdy(i,j,k) = KPV(i,j,k)*gradqy(i,j,k)
915 m_bates 1.3 XimY(m,i,j) = XimY(m,i,j)
916     & - drF(k)*hfacS(i,j,k,bi,bj)
917 m_bates 1.13 & *Kdqdy(i,j,k)*modesS(m,i,j,k,bi,bj)
918 m_bates 1.1 ENDDO
919     ENDDO
920     ENDDO
921     ENDDO
922    
923     C Calculate Xi for Y direction at the south face
924     DO k=1,Nr
925     DO j=1-Oly,sNy+Oly
926     DO i=1-Olx,sNx+Olx
927     Xiy(i,j,k) = zeroRL
928     ENDDO
929     ENDDO
930     ENDDO
931     DO k=1,Nr
932     DO j=1-Oly,sNy+Oly
933     DO i=1-Olx,sNx+Olx
934     DO m=1,GM_K3D_NModes
935     Xiy(i,j,k) = Xiy(i,j,k)
936 m_bates 1.4 & + maskS(i,j,k,bi,bj)*XimY(m,i,j)*modesS(m,i,j,k,bi,bj)
937 m_bates 1.1 ENDDO
938     ENDDO
939     ENDDO
940     ENDDO
941    
942 m_bates 1.11 C ENDIF (.NOT. GM_K3D_smooth)
943 m_bates 1.1 ENDIF
944    
945    
946 m_bates 1.15 C Calculate the renormalisation factor
947     DO j=1-Oly,sNy+Oly
948     DO i=1-Olx,sNx+Olx
949     uInt(i,j)=zeroRL
950     vInt(i,j)=zeroRL
951     KdqdyInt(i,j)=zeroRL
952     KdqdxInt(i,j)=zeroRL
953     uKdqdyInt(i,j)=zeroRL
954     vKdqdxInt(i,j)=zeroRL
955     uXiyInt(i,j)=zeroRL
956     vXixInt(i,j)=zeroRL
957 m_bates 1.16 Renorm(i,j)=oneRL
958     RenormU(i,j)=oneRL
959     RenormV(i,j)=oneRL
960 m_bates 1.15 ENDDO
961     ENDDO
962     DO k=1,Nr
963     DO j=1-Oly,sNy+Oly-1
964     DO i=1-Olx,sNx+Olx-1
965     centreX = op5*(uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))
966     centreY = op5*(Kdqdy(i,j,k) +Kdqdy(i,j+1,k) )
967     C For the numerator
968     uInt(i,j) = uInt(i,j)
969     & + centreX*hfacC(i,j,k,bi,bj)*drF(k)
970     KdqdyInt(i,j) = KdqdyInt(i,j)
971     & + centreY*hfacC(i,j,k,bi,bj)*drF(k)
972     uKdqdyInt(i,j) = uKdqdyInt(i,j)
973     & + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
974     C For the denominator
975     centreY = op5*(Xiy(i,j,k) + Xiy(i,j+1,k))
976     uXiyInt(i,j) = uXiyInt(i,j)
977     & + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
978    
979     centreX = op5*(Kdqdx(i,j,k) +Kdqdx(i+1,j,k))
980     centreY = op5*(vVel(i,j,k,bi,bj)+vVel(i,j+1,k,bi,bj) )
981     C For the numerator
982     vInt(i,j) = vInt(i,j)
983     & + centreY*hfacC(i,j,k,bi,bj)*drF(k)
984     KdqdxInt(i,j) = KdqdxInt(i,j)
985     & + CentreX*hfacC(i,j,k,bi,bj)*drF(k)
986     vKdqdxInt(i,j) = vKdqdxInt(i,j)
987     & + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
988     C For the denominator
989     centreX = op5*(Xix(i,j,k) + Xix(i+1,j,k))
990     vXixInt(i,j) = vXixInt(i,j)
991     & + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
992    
993     ENDDO
994     ENDDO
995     ENDDO
996    
997     DO j=1-Oly,sNy+Oly-1
998     DO i=1-Olx,sNx+Olx-1
999     IF (kLowC(i,j,bi,bj).GT.0) THEN
1000     numerator =
1001     & (uKdqdyInt(i,j)-uInt(i,j)*KdqdyInt(i,j)/R_low(i,j,bi,bj))
1002     & -(vKdqdxInt(i,j)-vInt(i,j)*KdqdxInt(i,j)/R_low(i,j,bi,bj))
1003     denominator = uXiyInt(i,j) - vXixInt(i,j)
1004     C We can have troubles with floating point exceptions if the denominator
1005     C of the renormalisation if the ocean is resting (e.g. intial conditions).
1006 m_bates 1.16 C So we make the renormalisation factor one if the denominator is very small
1007     C The renormalisation factor is supposed to correct the error in the extraction of
1008     C potential energy associated with the truncation of the expansion. Thus, we
1009     C enforce a minimum value for the renormalisation factor.
1010     C We also enforce a maximum renormalisation factor.
1011     IF (denominator.GT.small) THEN
1012 m_bates 1.15 Renorm(i,j) = ABS(numerator/denominator)
1013 m_bates 1.16 Renorm(i,j) = MAX(Renorm(i,j),GM_K3D_minRenorm)
1014     Renorm(i,j) = MIN(Renorm(i,j),GM_K3D_maxRenorm)
1015 m_bates 1.15 ENDIF
1016     ENDIF
1017     ENDDO
1018     ENDDO
1019     C Now put it back on to the velocity grids
1020     DO j=1-Oly+1,sNy+Oly-1
1021     DO i=1-Olx+1,sNx+Olx-1
1022     RenormU(i,j) = op5*(Renorm(i-1,j)+Renorm(i,j))
1023     RenormV(i,j) = op5*(Renorm(i,j-1)+Renorm(i,j))
1024     ENDDO
1025     ENDDO
1026    
1027 m_bates 1.1 C Calculate the eddy induced velocity in the X direction at the west face
1028     DO k=1,Nr
1029     DO j=1-Oly+1,sNy+Oly
1030     DO i=1-Olx+1,sNx+Olx
1031 m_bates 1.16 ustar(i,j,k) = -RenormU(i,j)*Xix(i,j,k)/coriU(i,j)
1032 m_bates 1.1 ENDDO
1033     ENDDO
1034     ENDDO
1035    
1036     C Calculate the eddy induced velocity in the Y direction at the south face
1037     DO k=1,Nr
1038     DO j=1-Oly+1,sNy+Oly
1039     DO i=1-Olx+1,sNx+Olx
1040 m_bates 1.16 vstar(i,j,k) = -RenormV(i,j)*Xiy(i,j,k)/coriV(i,j)
1041 m_bates 1.1 ENDDO
1042     ENDDO
1043     ENDDO
1044    
1045     C ======================================
1046     C Calculate the eddy induced overturning streamfunction
1047     C ======================================
1048     #ifdef GM_K3D_PASSIVE
1049     k=Nr
1050     DO j=1-Oly,sNy+Oly
1051     DO i=1-Olx,sNx+Olx
1052     psistar(i,j,Nr) = -hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
1053     ENDDO
1054     ENDDO
1055     DO k=Nr-1,1,-1
1056     DO j=1-Oly,sNy+Oly
1057     DO i=1-Olx,sNx+Olx
1058     psistar(i,j,k) = psistar(i,j,k+1)
1059     & - hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
1060     ENDDO
1061     ENDDO
1062     ENDDO
1063    
1064     #else
1065    
1066     IF (GM_AdvForm) THEN
1067     k=Nr
1068     DO j=1-Oly+1,sNy+1
1069     DO i=1-Olx+1,sNx+1
1070     GM_PsiX(i,j,k,bi,bj) = -hfacW(i,j,k,bi,bj)*drF(k)*ustar(i,j,k)
1071     GM_PsiY(i,j,k,bi,bj) = -hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
1072     ENDDO
1073     ENDDO
1074     DO k=Nr-1,1,-1
1075     DO j=1-Oly+1,sNy+1
1076     DO i=1-Olx+1,sNx+1
1077     GM_PsiX(i,j,k,bi,bj) = GM_PsiX(i,j,k+1,bi,bj)
1078     & - hfacW(i,j,k,bi,bj)*drF(k)*ustar(i,j,k)
1079     GM_PsiY(i,j,k,bi,bj) = GM_PsiY(i,j,k+1,bi,bj)
1080     & - hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
1081     ENDDO
1082     ENDDO
1083     ENDDO
1084    
1085     ENDIF
1086     #endif
1087    
1088     #ifdef ALLOW_DIAGNOSTICS
1089     C Diagnostics
1090     IF ( useDiagnostics ) THEN
1091     CALL DIAGNOSTICS_FILL(K3D, 'GM_K3D ',0,Nr,0,1,1,myThid)
1092 m_bates 1.14 CALL DIAGNOSTICS_FILL(KPV, 'GM_KPV ',0,Nr,0,1,1,myThid)
1093 m_bates 1.1 CALL DIAGNOSTICS_FILL(urms, 'GM_URMS ',0,Nr,0,1,1,myThid)
1094     CALL DIAGNOSTICS_FILL(Rdef, 'GM_RDEF ',0, 1,0,1,1,myThid)
1095 m_bates 1.8 CALL DIAGNOSTICS_FILL(Rurms, 'GM_RURMS',0, 1,0,1,1,myThid)
1096     CALL DIAGNOSTICS_FILL(RRhines,'GM_RRHNS',0, 1,0,1,1,myThid)
1097 m_bates 1.1 CALL DIAGNOSTICS_FILL(Rmix, 'GM_RMIX ',0, 1,0,1,1,myThid)
1098     CALL DIAGNOSTICS_FILL(supp, 'GM_SUPP ',0,Nr,0,1,1,myThid)
1099     CALL DIAGNOSTICS_FILL(Xix, 'GM_Xix ',0,Nr,0,1,1,myThid)
1100     CALL DIAGNOSTICS_FILL(Xiy, 'GM_Xiy ',0,Nr,0,1,1,myThid)
1101     CALL DIAGNOSTICS_FILL(cDopp, 'GM_C ',0, 1,0,1,1,myThid)
1102     CALL DIAGNOSTICS_FILL(Ubaro, 'GM_UBARO',0, 1,0,1,1,myThid)
1103     CALL DIAGNOSTICS_FILL(eady, 'GM_EADY ',0, 1,0,1,1,myThid)
1104     CALL DIAGNOSTICS_FILL(SlopeX, 'GM_Sx ',0,Nr,0,1,1,myThid)
1105     CALL DIAGNOSTICS_FILL(SlopeY, 'GM_Sy ',0,Nr,0,1,1,myThid)
1106     CALL DIAGNOSTICS_FILL(tfluxX, 'GM_TFLXX',0,Nr,0,1,1,myThid)
1107     CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,0,1,1,myThid)
1108     CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,0,1,1,myThid)
1109     CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,0,1,1,myThid)
1110 m_bates 1.13 CALL DIAGNOSTICS_FILL(Kdqdy, 'GM_Kdqdy',0,Nr,0,1,1,myThid)
1111     CALL DIAGNOSTICS_FILL(Kdqdx, 'GM_Kdqdx',0,Nr,0,1,1,myThid)
1112 m_bates 1.1 CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,0,1,1,myThid)
1113     CALL DIAGNOSTICS_FILL(ustar, 'GM_USTAR',0,Nr,0,1,1,myThid)
1114     CALL DIAGNOSTICS_FILL(vstar, 'GM_VSTAR',0,Nr,0,1,1,myThid)
1115     CALL DIAGNOSTICS_FILL(umc, 'GM_UMC ',0,Nr,0,1,1,myThid)
1116     CALL DIAGNOSTICS_FILL(ubar, 'GM_UBAR ',0,Nr,0,1,1,myThid)
1117 m_bates 1.4 CALL DIAGNOSTICS_FILL(modesC(1,:,:,:,bi,bj),
1118     & 'GM_MODEC',0,Nr,0,1,1,myThid)
1119 m_bates 1.10 CALL DIAGNOSTICS_FILL(M4loc, 'GM_M4 ',0,Nr,0,1,1,myThid)
1120     CALL DIAGNOSTICS_FILL(N2loc, 'GM_N2 ',0,Nr,0,1,1,myThid)
1121 m_bates 1.13 CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,0,1,1,myThid)
1122 m_bates 1.10 CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,0,1,1,myThid)
1123 m_bates 1.15 CALL DIAGNOSTICS_FILL(Renorm, 'GM_RENRM',0, 1,0,1,1,myThid)
1124 m_bates 1.10
1125 m_bates 1.1 ENDIF
1126     #endif
1127    
1128 m_bates 1.11 C For the Redi diffusivity, we set K3D to a constant if
1129 m_bates 1.14 C GM_K3D_constRedi=.TRUE.
1130     IF (GM_K3D_constRedi) THEN
1131 m_bates 1.11 DO k=1,Nr
1132     DO j=1-Oly,sNy+Oly
1133     DO i=1-Olx,sNx+Olx
1134     K3D(i,j,k,bi,bj) = GM_K3D_constK
1135     ENDDO
1136     ENDDO
1137     ENDDO
1138     ENDIF
1139    
1140 m_bates 1.14 #ifdef ALLOW_DIAGNOSTICS
1141     IF ( useDiagnostics )
1142     & CALL DIAGNOSTICS_FILL(K3D, 'GM_K3D_T',0,Nr,0,1,1,myThid)
1143     #endif
1144    
1145 m_bates 1.1 #endif /* GM_K3D */
1146     RETURN
1147     END
  ViewVC Help
Powered by ViewVC 1.1.22