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