ROMS
Loading...
Searching...
No Matches
my25_corstep.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE my25_corstep_mod
3#if defined NONLINEAR && defined MY25_MIXING && defined SOLVE3D
4!
5!git $Id$
6!=======================================================================
7! Copyright (c) 2002-2025 The ROMS Group !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md Hernan G. Arango !
10!========================================== Alexander F. Shchepetkin ===
11! !
12! This routine perfoms the corrector step for turbulent kinetic !
13! energy and length scale prognostic variables, tke and gls. !
14! !
15! References: !
16! !
17! Mellor, G.L. and T. Yamada, 1982: Development of turbulence !
18! closure model for geophysical fluid problems, Rev. Geophys. !
19! Space Phys., 20, 851-875. !
20! !
21! Galperin, B., L.H. Kantha, S. Hassid, and A.Rosati, 1988: A !
22! quasi-equilibrium turbulent energy model for geophysical !
23! flows, J. Atmos. Sci., 45, 55-62. !
24! !
25!=======================================================================
26!
27 implicit none
28!
29 PRIVATE
30 PUBLIC :: my25_corstep
31!
32 CONTAINS
33!
34!***********************************************************************
35 SUBROUTINE my25_corstep (ng, tile)
36!***********************************************************************
37!
38 USE mod_param
39 USE mod_forces
40 USE mod_grid
41 USE mod_mixing
42 USE mod_ocean
43 USE mod_stepping
44!
45! Imported variable declarations.
46!
47 integer, intent(in) :: ng, tile
48!
49! Local variable declarations.
50!
51 character (len=*), parameter :: myfile = &
52 & __FILE__
53!
54# include "tile.h"
55!
56# ifdef PROFILE
57 CALL wclock_on (ng, inlm, 20, __line__, myfile)
58# endif
59 CALL my25_corstep_tile (ng, tile, &
60 & lbi, ubi, lbj, ubj, &
61 & imins, imaxs, jmins, jmaxs, &
62 & nstp(ng), nnew(ng), &
63# ifdef MASKING
64 & grid(ng) % umask, &
65 & grid(ng) % vmask, &
66# endif
67 & grid(ng) % Huon, &
68 & grid(ng) % Hvom, &
69 & grid(ng) % Hz, &
70 & grid(ng) % pm, &
71 & grid(ng) % pn, &
72 & grid(ng) % z_r, &
73 & grid(ng) % z_w, &
74 & ocean(ng) % u, &
75 & ocean(ng) % v, &
76 & ocean(ng) % W, &
77 & forces(ng) % bustr, &
78 & forces(ng) % bvstr, &
79 & forces(ng) % sustr, &
80 & forces(ng) % svstr, &
81 & mixing(ng) % bvf, &
82 & mixing(ng) % Akt, &
83 & mixing(ng) % Akv, &
84 & mixing(ng) % Akk, &
85 & mixing(ng) % Lscale, &
86 & mixing(ng) % gls, &
87 & mixing(ng) % tke)
88# ifdef PROFILE
89 CALL wclock_off (ng, inlm, 20, __line__, myfile)
90# endif
91!
92 RETURN
93 END SUBROUTINE my25_corstep
94!
95!***********************************************************************
96 SUBROUTINE my25_corstep_tile (ng, tile, &
97 & LBi, UBi, LBj, UBj, &
98 & IminS, ImaxS, JminS, JmaxS, &
99 & nstp, nnew, &
100# ifdef MASKING
101 & umask, vmask, &
102# endif
103 & Huon, Hvom, Hz, pm, pn, z_r, z_w, &
104 & u, v, W, &
105 & bustr, bvstr, sustr, svstr, &
106 & bvf, Akt, Akv, &
107 & Akk, Lscale, gls, tke)
108!***********************************************************************
109!
110 USE mod_param
111 USE mod_ncparam
112 USE mod_scalars
113!
114 USE exchange_3d_mod, ONLY : exchange_w3d_tile
115# ifdef DISTRIBUTE
116 USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
117# endif
118 USE tkebc_mod, ONLY : tkebc_tile
119!
120! Imported variable declarations.
121!
122 integer, intent(in) :: ng, tile
123 integer, intent(in) :: LBi, UBi, LBj, UBj
124 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
125 integer, intent(in) :: nstp, nnew
126!
127# ifdef ASSUMED_SHAPE
128# ifdef MASKING
129 real(r8), intent(in) :: umask(LBi:,LBj:)
130 real(r8), intent(in) :: vmask(LBi:,LBj:)
131# endif
132 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
133 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
134 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
135 real(r8), intent(in) :: pm(LBi:,LBj:)
136 real(r8), intent(in) :: pn(LBi:,LBj:)
137 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
138 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
139 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
140 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
141 real(r8), intent(in) :: W(LBi:,LBj:,0:)
142 real(r8), intent(in) :: bustr(LBi:,LBj:)
143 real(r8), intent(in) :: bvstr(LBi:,LBj:)
144 real(r8), intent(in) :: sustr(LBi:,LBj:)
145 real(r8), intent(in) :: svstr(LBi:,LBj:)
146 real(r8), intent(in) :: bvf(LBi:,LBj:,0:)
147
148 real(r8), intent(inout) :: Akt(LBi:,LBj:,0:,:)
149 real(r8), intent(inout) :: Akv(LBi:,LBj:,0:)
150 real(r8), intent(inout) :: Akk(LBi:,LBj:,0:)
151 real(r8), intent(inout) :: Lscale(LBi:,LBj:,0:)
152 real(r8), intent(inout) :: gls(LBi:,LBj:,0:,:)
153 real(r8), intent(inout) :: tke(LBi:,LBj:,0:,:)
154# else
155# ifdef MASKING
156 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
157 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
158# endif
159 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
160 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
161 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
162 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
163 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
164 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
165 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
166 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
167 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
168 real(r8), intent(in) :: W(LBi:UBi,LBj:UBj,0:N(ng))
169 real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj)
170 real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
171 real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj)
172 real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj)
173 real(r8), intent(in) :: bvf(LBi:UBi,LBj:UBj,0:N(ng))
174
175 real(r8), intent(inout) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
176 real(r8), intent(inout) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
177 real(r8), intent(inout) :: Akk(LBi:UBi,LBj:UBj,0:N(ng))
178 real(r8), intent(inout) :: Lscale(LBi:UBi,LBj:UBj,0:N(ng))
179 real(r8), intent(inout) :: gls(LBi:UBi,LBj:UBj,0:N(ng),3)
180 real(r8), intent(inout) :: tke(LBi:UBi,LBj:UBj,0:N(ng),3)
181# endif
182!
183! Local variable declarations.
184!
185 integer :: i, itrc, j, k
186
187 real(r8), parameter :: Gadv = 1.0_r8/3.0_r8
188 real(r8), parameter :: eps = 1.0e-10_r8
189
190 real(r8) :: Gh, Ls_unlmt, Ls_lmt, Qprod, Qdiss, Sh, Sm, Wscale
191 real(r8) :: cff, cff1, cff2, cff3, ql, strat2
192
193 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BCK
194 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BCP
195 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
196 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCK
197 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCP
198 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dU
199 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dV
200
201 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: shear2
202 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: buoy2
203
204 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FEK
205 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FEP
206 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FXK
207 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FXP
208 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curvK
209 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curvP
210 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradK
211 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradP
212
213# include "set_bounds.h"
214!
215!-----------------------------------------------------------------------
216! Compute vertical velocity shear at W-points.
217!-----------------------------------------------------------------------
218!
219# ifdef RI_SPLINES
220 DO j=jstrm1,jendp1
221 DO i=istrm1,iendp1
222 cf(i,0)=0.0_r8
223 du(i,0)=0.0_r8
224 dv(i,0)=0.0_r8
225 END DO
226 DO k=1,n(ng)-1
227 DO i=istrm1,iendp1
228 cff=1.0_r8/(2.0_r8*hz(i,j,k+1)+ &
229 & hz(i,j,k)*(2.0_r8-cf(i,k-1)))
230 cf(i,k)=cff*hz(i,j,k+1)
231 du(i,k)=cff*(3.0_r8*(u(i ,j,k+1,nstp)-u(i, j,k,nstp)+ &
232 & u(i+1,j,k+1,nstp)-u(i+1,j,k,nstp))- &
233 & hz(i,j,k)*du(i,k-1))
234 dv(i,k)=cff*(3.0_r8*(v(i,j ,k+1,nstp)-v(i,j ,k,nstp)+ &
235 & v(i,j+1,k+1,nstp)-v(i,j+1,k,nstp))- &
236 & hz(i,j,k)*dv(i,k-1))
237 END DO
238 END DO
239 DO i=istrm1,iendp1
240 du(i,n(ng))=0.0_r8
241 dv(i,n(ng))=0.0_r8
242 END DO
243 DO k=n(ng)-1,1,-1
244 DO i=istrm1,iendp1
245 du(i,k)=du(i,k)-cf(i,k)*du(i,k+1)
246 dv(i,k)=dv(i,k)-cf(i,k)*dv(i,k+1)
247 END DO
248 END DO
249 DO k=1,n(ng)-1
250 DO i=istrm1,iendp1
251 shear2(i,j,k)=du(i,k)*du(i,k)+dv(i,k)*dv(i,k)
252 END DO
253 END DO
254 END DO
255# else
256 DO k=1,n(ng)-1
257 DO j=jstrm1,jendp1
258 DO i=istrm1,iendp1
259 cff=0.5_r8/(z_r(i,j,k+1)-z_r(i,j,k))
260 shear2(i,j,k)=(cff*(u(i ,j,k+1,nstp)-u(i ,j,k,nstp)+ &
261 & u(i+1,j,k+1,nstp)-u(i+1,j,k,nstp)))**2+ &
262 & (cff*(v(i,j ,k+1,nstp)-v(i,j ,k,nstp)+ &
263 & v(i,j+1,k+1,nstp)-v(i,j+1,k,nstp)))**2
264 END DO
265 END DO
266 END DO
267# endif
268!
269! Load Brunt-Vaisala frequency.
270!
271 DO k=1,n(ng)-1
272 DO j=jstr-1,jend+1
273 DO i=istr-1,iend+1
274 buoy2(i,j,k)=bvf(i,j,k)
275 END DO
276 END DO
277 END DO
278# ifdef N2S2_HORAVG
279!
280!-----------------------------------------------------------------------
281! Smooth horizontally buoyancy and shear. Use buoy2(:,:,0) and
282! shear2(:,:,0) as scratch utility array.
283!-----------------------------------------------------------------------
284!
285 DO k=1,n(ng)-1
286 IF (domain(ng)%Western_Edge(tile)) THEN
287 DO j=max(1,jstr-1),min(jend+1,mm(ng))
288 shear2(istr-1,j,k)=shear2(istr,j,k)
289 END DO
290 END IF
291 IF (domain(ng)%Eastern_Edge(tile)) THEN
292 DO j=max(1,jstr-1),min(jend+1,mm(ng))
293 shear2(iend+1,j,k)=shear2(iend,j,k)
294 END DO
295 END IF
296 IF (domain(ng)%Southern_Edge(tile)) THEN
297 DO i=max(1,istr-1),min(iend+1,lm(ng))
298 shear2(i,jstr-1,k)=shear2(i,jstr,k)
299 END DO
300 END IF
301 IF (domain(ng)%Northern_Edge(tile)) THEN
302 DO i=max(1,istr-1),min(iend+1,lm(ng))
303 shear2(i,jend+1,k)=shear2(i,jend,k)
304 END DO
305 END IF
306 IF (domain(ng)%SouthWest_Corner(tile)) THEN
307 shear2(istr-1,jstr-1,k)=shear2(istr,jstr,k)
308 END IF
309 IF (domain(ng)%NorthWest_Corner(tile)) THEN
310 shear2(istr-1,jend+1,k)=shear2(istr,jend,k)
311 END IF
312 IF (domain(ng)%SouthEast_Corner(tile)) THEN
313 shear2(iend+1,jstr-1,k)=shear2(iend,jstr,k)
314 END IF
315 IF (domain(ng)%NorthEast_Corner(tile)) THEN
316 shear2(iend+1,jend+1,k)=shear2(iend,jend,k)
317 END IF
318!
319! Average horizontally.
320!
321 DO j=jstr-1,jend
322 DO i=istr-1,iend
323 buoy2(i,j,0)=0.25_r8*(buoy2(i,j ,k)+buoy2(i+1,j ,k)+ &
324 & buoy2(i,j+1,k)+buoy2(i+1,j+1,k))
325 shear2(i,j,0)=0.25_r8*(shear2(i,j ,k)+shear2(i+1,j ,k)+ &
326 & shear2(i,j+1,k)+shear2(i+1,j+1,k))
327 END DO
328 END DO
329 DO j=jstr,jend
330 DO i=istr,iend
331 buoy2(i,j,k)=0.25_r8*(buoy2(i,j ,0)+buoy2(i-1,j ,0)+ &
332 & buoy2(i,j-1,0)+buoy2(i-1,j-1,0))
333 shear2(i,j,k)=0.25_r8*(shear2(i,j ,0)+shear2(i-1,j ,0)+ &
334 & shear2(i,j-1,0)+shear2(i-1,j-1,0))
335 END DO
336 END DO
337 END DO
338# endif
339!
340!-----------------------------------------------------------------------
341! Time-step advective terms.
342!-----------------------------------------------------------------------
343!
344! At entry, it is assumed that the turbulent kinetic energy fields
345! "tke" and "gls", at time level "nnew", are set to its values at
346! time level "nstp" multiplied by the grid box thicknesses Hz
347! (from old time step and at W-points).
348!
349 DO k=1,n(ng)-1
350# ifdef K_C2ADVECTION
351!
352! Second-order, centered differences advection.
353!
354 DO j=jstr,jend
355 DO i=istr,iend+1
356 cff=0.25_r8*(huon(i,j,k)+huon(i,j,k+1))
357 fxk(i,j)=cff*(tke(i,j,k,3)+tke(i-1,j,k,3))
358 fxp(i,j)=cff*(gls(i,j,k,3)+gls(i-1,j,k,3))
359 END DO
360 END DO
361 DO j=jstr,jend+1
362 DO i=istr,iend
363 cff=0.25_r8*(hvom(i,j,k)+hvom(i,j,k+1))
364 fek(i,j)=cff*(tke(i,j,k,3)+tke(i,j-1,k,3))
365 fep(i,j)=cff*(gls(i,j,k,3)+gls(i,j-1,k,3))
366 END DO
367 END DO
368# else
369 DO j=jstr,jend
370 DO i=istrm1,iendp2
371 gradk(i,j)=(tke(i,j,k,3)-tke(i-1,j,k,3))
372# ifdef MASKING
373 gradk(i,j)=gradk(i,j)*umask(i,j)
374# endif
375 gradp(i,j)=(gls(i,j,k,3)-gls(i-1,j,k,3))
376# ifdef MASKING
377 gradp(i,j)=gradp(i,j)*umask(i,j)
378# endif
379 END DO
380 END DO
381 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
382 IF (domain(ng)%Western_Edge(tile)) THEN
383 DO j=jstr,jend
384 gradk(istr-1,j)=gradk(istr,j)
385 gradp(istr-1,j)=gradp(istr,j)
386 END DO
387 END IF
388 END IF
389 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
390 IF (domain(ng)%Eastern_Edge(tile)) THEN
391 DO j=jstr,jend
392 gradk(iend+2,j)=gradk(iend+1,j)
393 gradp(iend+2,j)=gradp(iend+1,j)
394 END DO
395 END IF
396 END IF
397# ifdef K_C4ADVECTION
398!
399! Fourth-order, centered differences advection.
400!
401 cff1=1.0_r8/6.0_r8
402 DO j=jstr,jend
403 DO i=istr,iend+1
404 cff=0.5_r8*(huon(i,j,k)+huon(i,j,k+1))
405 fxk(i,j)=cff*0.5_r8*(tke(i-1,j,k,3)+tke(i,j,k,3)- &
406 & cff1*(gradk(i+1,j)-gradk(i-1,j)))
407 fxp(i,j)=cff*0.5_r8*(gls(i-1,j,k,3)+gls(i,j,k,3)- &
408 & cff1*(gradp(i+1,j)-gradp(i-1,j)))
409 END DO
410 END DO
411# else
412!
413! Third-order, upstream bias advection with velocity dependent
414! hyperdiffusion.
415!
416 DO j=jstr,jend
417 DO i=istr-1,iend+1
418 curvk(i,j)=gradk(i+1,j)-gradk(i,j)
419 curvp(i,j)=gradp(i+1,j)-gradp(i,j)
420 END DO
421 END DO
422 DO j=jstr,jend
423 DO i=istr,iend+1
424 cff=0.5_r8*(huon(i,j,k)+huon(i,j,k+1))
425 IF (cff.gt.0.0_r8) THEN
426 cff1=curvk(i-1,j)
427 cff2=curvp(i-1,j)
428 ELSE
429 cff1=curvk(i,j)
430 cff2=curvp(i,j)
431 END IF
432 fxk(i,j)=cff*0.5_r8*(tke(i-1,j,k,3)+tke(i,j,k,3)- &
433 & gadv*cff1)
434 fxp(i,j)=cff*0.5_r8*(gls(i-1,j,k,3)+gls(i,j,k,3)- &
435 & gadv*cff2)
436 END DO
437 END DO
438# endif
439 DO j=jstrm1,jendp2
440 DO i=istr,iend
441 gradk(i,j)=(tke(i,j,k,3)-tke(i,j-1,k,3))
442# ifdef MASKING
443 gradk(i,j)=gradk(i,j)*vmask(i,j)
444# endif
445 gradp(i,j)=(gls(i,j,k,3)-gls(i,j-1,k,3))
446# ifdef MASKING
447 gradp(i,j)=gradp(i,j)*vmask(i,j)
448# endif
449 END DO
450 END DO
451 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
452 IF (domain(ng)%Southern_Edge(tile)) THEN
453 DO i=istr,iend
454 gradk(i,jstr-1)=gradk(i,jstr)
455 gradp(i,jstr-1)=gradp(i,jstr)
456 END DO
457 END IF
458 END IF
459 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
460 IF (domain(ng)%Northern_Edge(tile)) THEN
461 DO i=istr,iend
462 gradk(i,jend+2)=gradk(i,jend+1)
463 gradp(i,jend+2)=gradp(i,jend+1)
464 END DO
465 END IF
466 END IF
467# ifdef K_C4ADVECTION
468 cff1=1.0_r8/6.0_r8
469 DO j=jstr,jend+1
470 DO i=istr,iend
471 cff=0.5_r8*(hvom(i,j,k)+hvom(i,j,k+1))
472 fek(i,j)=cff*0.5_r8*(tke(i,j-1,k,3)+tke(i,j,k,3)- &
473 & cff1*(gradk(i,j+1)-gradk(i,j-1)))
474 fep(i,j)=cff*0.5_r8*(gls(i,j-1,k,3)+gls(i,j,k,3)- &
475 & cff1*(gradp(i,j+1)-gradp(i,j-1)))
476 END DO
477 END DO
478# else
479 DO j=jstr-1,jend+1
480 DO i=istr,iend
481 curvk(i,j)=gradk(i,j+1)-gradk(i,j)
482 curvp(i,j)=gradp(i,j+1)-gradp(i,j)
483 END DO
484 END DO
485 DO j=jstr,jend+1
486 DO i=istr,iend
487 cff=0.5_r8*(hvom(i,j,k)+hvom(i,j,k+1))
488 IF (cff.gt.0.0_r8) THEN
489 cff1=curvk(i,j-1)
490 cff2=curvp(i,j-1)
491 ELSE
492 cff1=curvk(i,j)
493 cff2=curvp(i,j)
494 END IF
495 fek(i,j)=cff*0.5_r8*(tke(i,j-1,k,3)+tke(i,j,k,3)- &
496 & gadv*cff1)
497 fep(i,j)=cff*0.5_r8*(gls(i,j-1,k,3)+gls(i,j,k,3)- &
498 & gadv*cff2)
499 END DO
500 END DO
501# endif
502# endif
503!
504! Time-step horizontal advection.
505!
506 DO j=jstr,jend
507 DO i=istr,iend
508 cff=dt(ng)*pm(i,j)*pn(i,j)
509 tke(i,j,k,nnew)=tke(i,j,k,nnew)- &
510 & cff*(fxk(i+1,j)-fxk(i,j)+ &
511 & fek(i,j+1)-fek(i,j))
512 gls(i,j,k,nnew)=gls(i,j,k,nnew)- &
513 & cff*(fxp(i+1,j)-fxp(i,j)+ &
514 & fep(i,j+1)-fep(i,j))
515 END DO
516 END DO
517 END DO
518!
519! Compute vertical advection.
520!
521 DO j=jstr,jend
522# ifdef K_C2ADVECTION
523 DO k=1,n(ng)
524 DO i=istr,iend
525 cff=0.25_r8*(w(i,j,k)+w(i,j,k-1))
526 fck(i,k)=cff*(tke(i,j,k,3)+tke(i,j,k-1,3))
527 fcp(i,k)=cff*(gls(i,j,k,3)+gls(i,j,k-1,3))
528 END DO
529 END DO
530# else
531 cff1=7.0_r8/12.0_r8
532 cff2=1.0_r8/12.0_r8
533 DO k=2,n(ng)-1
534 DO i=istr,iend
535 cff=0.5*(w(i,j,k)+w(i,j,k-1))
536 fck(i,k)=cff*(cff1*(tke(i,j,k-1,3)+ &
537 & tke(i,j,k ,3))- &
538 & cff2*(tke(i,j,k-2,3)+ &
539 & tke(i,j,k+1,3)))
540 fcp(i,k)=cff*(cff1*(gls(i,j,k-1,3)+ &
541 & gls(i,j,k ,3))- &
542 & cff2*(gls(i,j,k-2,3)+ &
543 & gls(i,j,k+1,3)))
544 END DO
545 END DO
546 cff1=1.0_r8/3.0_r8
547 cff2=5.0_r8/6.0_r8
548 cff3=1.0_r8/6.0_r8
549 DO i=istr,iend
550 cff=0.5_r8*(w(i,j,0)+w(i,j,1))
551 fck(i,1)=cff*(cff1*tke(i,j,0,3)+ &
552 & cff2*tke(i,j,1,3)- &
553 & cff3*tke(i,j,2,3))
554 fcp(i,1)=cff*(cff1*gls(i,j,0,3)+ &
555 & cff2*gls(i,j,1,3)- &
556 & cff3*gls(i,j,2,3))
557 cff=0.5_r8*(w(i,j,n(ng))+w(i,j,n(ng)-1))
558 fck(i,n(ng))=cff*(cff1*tke(i,j,n(ng) ,3)+ &
559 & cff2*tke(i,j,n(ng)-1,3)- &
560 & cff3*tke(i,j,n(ng)-2,3))
561 fcp(i,n(ng))=cff*(cff1*gls(i,j,n(ng) ,3)+ &
562 & cff2*gls(i,j,n(ng)-1,3)- &
563 & cff3*gls(i,j,n(ng)-2,3))
564 END DO
565# endif
566!
567! Time-step vertical advection term.
568!
569 DO k=1,n(ng)-1
570 DO i=istr,iend
571 cff=dt(ng)*pm(i,j)*pn(i,j)
572 tke(i,j,k,nnew)=tke(i,j,k,nnew)- &
573 & cff*(fck(i,k+1)-fck(i,k))
574 gls(i,j,k,nnew)=gls(i,j,k,nnew)- &
575 & cff*(fcp(i,k+1)-fcp(i,k))
576 END DO
577 END DO
578!
579!----------------------------------------------------------------------
580! Compute vertical mixing, turbulent production and turbulent
581! dissipation terms.
582!----------------------------------------------------------------------
583!
584! Set term for vertical mixing of turbulent fields.
585!
586 cff=-0.5_r8*dt(ng)
587 DO k=1,n(ng)
588 DO i=istr,iend
589 fck(i,k)=cff*(akk(i,j,k)+akk(i,j,k-1))/hz(i,j,k)
590 cf(i,k)=0.0_r8
591 END DO
592 END DO
593!
594! Compute production and dissipation terms.
595!
596 cff3=my_e2/(vonkar*vonkar)
597 DO k=1,n(ng)-1
598 DO i=istr,iend
599!
600! Compute shear and bouyant production of turbulent energy (m3/s3)
601! at W-points (ignore small negative values of buoyancy).
602!
603 IF ((buoy2(i,j,k).gt.-5.0e-5_r8).and. &
604 & (buoy2(i,j,k).lt.0.0_r8)) THEN
605 strat2=0.0_r8
606 ELSE
607 strat2=buoy2(i,j,k)
608 END IF
609 qprod=shear2(i,j,k)*(akv(i,j,k)-akv_bak(ng))- &
610 & strat2*(akt(i,j,k,itemp)-akt_bak(itemp,ng))
611!
612! Recalculate old time-step unlimited length scale.
613!
614 ls_unlmt=max(eps, &
615 & gls(i,j,k,nstp)/(max(tke(i,j,k,nstp),eps)))
616!
617! Time-step production term.
618!
619 cff1=0.5_r8*(hz(i,j,k)+hz(i,j,k+1))
620 tke(i,j,k,nnew)=tke(i,j,k,nnew)+ &
621 & dt(ng)*cff1*qprod*2.0_r8
622 gls(i,j,k,nnew)=gls(i,j,k,nnew)+ &
623 & dt(ng)*cff1*qprod*my_e1*ls_unlmt
624!
625! Compute dissipation of turbulent energy (m3/s3). Add in vertical
626! mixing term.
627!
628 qdiss=dt(ng)*sqrt(tke(i,j,k,nstp))/(my_b1*ls_unlmt)
629 cff=ls_unlmt*(1.0_r8/(z_w(i,j,n(ng))-z_w(i,j,k))+ &
630 & 1.0_r8/(z_w(i,j,k)-z_w(i,j,0)))
631 wscale=1.0_r8+cff3*cff*cff
632 bck(i,k)=cff1*(1.0_r8+2.0_r8*qdiss)-fck(i,k)-fck(i,k+1)
633 bcp(i,k)=cff1*(1.0_r8+wscale*qdiss)-fck(i,k)-fck(i,k+1)
634 END DO
635 END DO
636!
637!-----------------------------------------------------------------------
638! Time-step dissipation and vertical diffusion terms implicitly.
639!-----------------------------------------------------------------------
640!
641! Set surface and bottom boundary conditions.
642!
643 DO i=istr,iend
644 tke(i,j,n(ng),nnew)=my_b1p2o3*0.5_r8* &
645 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
646 & (svstr(i,j)+svstr(i,j+1))**2)
647 gls(i,j,n(ng),nnew)=0.0_r8
648 tke(i,j,0,nnew)=my_b1p2o3*0.5_r8* &
649 & sqrt((bustr(i,j)+bustr(i+1,j))**2+ &
650 & (bvstr(i,j)+bvstr(i,j+1))**2)
651 gls(i,j,0,nnew)=0.0_r8
652 END DO
653!
654! Solve tri-diagonal system for "tke".
655!
656 DO i=istr,iend
657 cff=1.0_r8/bck(i,n(ng)-1)
658 cf(i,n(ng)-1)=cff*fck(i,n(ng)-1)
659 tke(i,j,n(ng)-1,nnew)=cff*(tke(i,j,n(ng)-1,nnew)- &
660 & fck(i,n(ng))*tke(i,j,n(ng),nnew))
661 END DO
662 DO k=n(ng)-2,1,-1
663 DO i=istr,iend
664 cff=1.0_r8/(bck(i,k)-cf(i,k+1)*fck(i,k+1))
665 cf(i,k)=cff*fck(i,k)
666 tke(i,j,k,nnew)=cff*(tke(i,j,k,nnew)- &
667 & fck(i,k+1)*tke(i,j,k+1,nnew))
668 END DO
669 END DO
670 DO k=1,n(ng)-1
671 DO i=istr,iend
672 tke(i,j,k,nnew)=tke(i,j,k,nnew)-cf(i,k)*tke(i,j,k-1,nnew)
673 END DO
674 END DO
675!
676! Solve tri-diagonal system for "gls".
677!
678 DO i=istr,iend
679 cff=1.0_r8/bcp(i,n(ng)-1)
680 cf(i,n(ng)-1)=cff*fck(i,n(ng)-1)
681 gls(i,j,n(ng)-1,nnew)=cff*(gls(i,j,n(ng)-1,nnew)- &
682 & fck(i,n(ng))*gls(i,j,n(ng),nnew))
683 END DO
684 DO k=n(ng)-2,1,-1
685 DO i=istr,iend
686 cff=1.0_r8/(bcp(i,k)-cf(i,k+1)*fck(i,k+1))
687 cf(i,k)=cff*fck(i,k)
688 gls(i,j,k,nnew)=cff*(gls(i,j,k,nnew)- &
689 & fck(i,k+1)*gls(i,j,k+1,nnew))
690 END DO
691 END DO
692 DO k=1,n(ng)-1,+1
693 DO i=istr,iend
694 gls(i,j,k,nnew)=gls(i,j,k,nnew)-cf(i,k)*gls(i,j,k-1,nnew)
695 END DO
696 END DO
697!
698!---------------------------------------------------------------------
699! Compute vertical mixing coefficients (m2/s).
700!---------------------------------------------------------------------
701!
702 DO k=1,n(ng)-1
703 DO i=istr,iend
704!
705! Compute turbulent length scale (m). The length scale is only
706! limited in the K-related calculations and not in QL production,
707! dissipation, wall-proximity, etc.
708!
709 tke(i,j,k,nnew)=max(tke(i,j,k,nnew),my_qmin)
710 gls(i,j,k,nnew)=max(gls(i,j,k,nnew),my_qmin)
711 ls_unlmt=gls(i,j,k,nnew)/tke(i,j,k,nnew)
712 ls_lmt=min(ls_unlmt, &
713 & my_lmax*sqrt(tke(i,j,k,nnew)/ &
714 & (max(0.0_r8,buoy2(i,j,k))+eps)))
715!
716! Compute Galperin et al. (1988) nondimensional stability function,
717! Gh. Then, compute nondimensional stability functions for tracers
718! (Sh) and momentum (Sm). The limit on length scale, sets the lower
719! limit on Gh. Use Kantha and Clayton or Galperin et al. expression
720! for Sm.
721!
722 gh=min(my_gh0,-buoy2(i,j,k)*ls_lmt*ls_lmt/ &
723 & tke(i,j,k,nnew))
724 cff=1.0_r8-my_sh2*gh
725 sh=my_sh1/cff
726# ifdef KANTHA_CLAYSON
727 sm=(my_b1pm1o3+sh*gh*my_sm4)/(1.0_r8-my_sm2*gh)
728# else
729 sm=(my_sm3+sh*gh*my_sm4)/(1.0_r8-my_sm2*gh)
730# endif
731!
732! Compute vertical mixing (m2/s) coefficients of momentum and
733! tracers. Average ql over the two timesteps rather than using
734! the new Lscale and just averaging tke.
735!
736 ql=0.5_r8*(ls_lmt*sqrt(tke(i,j,k,nnew))+ &
737 & lscale(i,j,k)*sqrt(tke(i,j,k,nstp)))
738 akv(i,j,k)=akv_bak(ng)+ql*sm
739 DO itrc=1,nat
740 akt(i,j,k,itrc)=akt_bak(itrc,ng)+ql*sh
741 END DO
742!
743! Compute vertical mixing (m2/s) coefficient of turbulent kinetic
744! energy. Use original formulation (Mellor and Yamada 1982;
745! Blumberg 1991; Kantha and Clayson 1994).
746!
747 akk(i,j,k)=akk_bak(ng)+ql*my_sq
748!
749! Save limited length scale.
750!
751 lscale(i,j,k)=ls_lmt
752
753# if defined LIMIT_VDIFF || defined LIMIT_VVISC
754!
755! Limit vertical mixing coefficients with the upper threshold value.
756! This is an engineering fix but it can be based on the fact that
757! vertical mixing in the ocean from indirect observations are not
758! higher than the threshold value.
759!
760# ifdef LIMIT_VDIFF
761 DO itrc=1,nat
762 akt(i,j,k,itrc)=min(akt_limit(itrc,ng), akt(i,j,k,itrc))
763 END DO
764# endif
765# ifdef LIMIT_VVISC
766 akv(i,j,k)=min(akv_limit(ng), akv(i,j,k))
767# endif
768# endif
769 END DO
770 END DO
771 END DO
772!
773!-----------------------------------------------------------------------
774! Set lateral boundary conditions.
775!-----------------------------------------------------------------------
776!
777 DO k=0,n(ng)
778 IF (domain(ng)%Western_Edge(tile)) THEN
779 DO j=jstr,jend
780 DO itrc=1,nat
781 akt(istr-1,j,k,itrc)=akt(istr,j,k,itrc)
782 END DO
783 akv(istr-1,j,k)=akv(istr,j,k)
784 END DO
785 END IF
786 IF (domain(ng)%Eastern_Edge(tile)) THEN
787 DO j=jstr,jend
788 DO itrc=1,nat
789 akt(iend-1,j,k,itrc)=akt(iend,j,k,itrc)
790 END DO
791 akv(iend-1,j,k)=akv(iend,j,k)
792 END DO
793
794 END IF
795 IF (domain(ng)%Southern_Edge(tile)) THEN
796 DO i=istr,iend
797 DO itrc=1,nat
798 akt(i,jstr-1,k,itrc)=akt(i,jstr,k,itrc)
799 END DO
800 akv(i,jstr-1,k)=akv(i,jstr,k)
801 END DO
802 END IF
803 IF (domain(ng)%Northern_Edge(tile)) THEN
804 DO i=istr,iend
805 DO itrc=1,nat
806 akt(i,jend+1,k,itrc)=akt(i,jend,k,itrc)
807 END DO
808 akv(i,jend+1,k)=akv(i,jend,k)
809 END DO
810 END IF
811 IF (domain(ng)%SouthWest_Corner(tile)) THEN
812 DO itrc=1,nat
813 akt(istr-1,jstr-1,k,itrc)=0.5_r8* &
814 & (akt(istr ,jstr-1,k,itrc)+ &
815 & akt(istr-1,jstr ,k,itrc))
816 END DO
817 akv(istr-1,jstr-1,k)=0.5_r8* &
818 & (akv(istr ,jstr-1,k)+ &
819 & akv(istr-1,jstr ,k))
820 END IF
821 IF (domain(ng)%SouthEast_Corner(tile)) THEN
822 DO itrc=1,nat
823 akt(iend+1,jstr-1,k,itrc)=0.5_r8* &
824 & (akt(iend ,jstr-1,k,itrc)+ &
825 & akt(iend+1,jstr ,k,itrc))
826 END DO
827 akv(iend+1,jstr-1,k)=0.5_r8* &
828 & (akv(iend ,jstr-1,k)+ &
829 & akv(iend+1,jstr ,k))
830 END IF
831 IF (domain(ng)%NorthWest_Corner(tile)) THEN
832 DO itrc=1,nat
833 akt(istr-1,jend+1,k,itrc)=0.5_r8* &
834 & (akt(istr ,jend+1,k,itrc)+ &
835 & akt(istr-1,jend ,k,itrc))
836 END DO
837 akv(istr-1,jend+1,k)=0.5_r8* &
838 & (akv(istr ,jend+1,k)+ &
839 & akv(istr-1,jend ,k))
840 END IF
841 IF (domain(ng)%NorthEast_Corner(tile)) THEN
842 DO itrc=1,nat
843 akt(iend+1,jend+1,k,itrc)=0.5_r8* &
844 & (akt(iend ,jend+1,k,itrc)+ &
845 & akt(iend+1,jend ,k,itrc))
846 END DO
847 akv(iend+1,jend+1,k)=0.5_r8* &
848 & (akv(iend ,jend+1,k)+ &
849 & akv(iend+1,jend ,k))
850 END IF
851 END DO
852!
853 CALL tkebc_tile (ng, tile, &
854 & lbi, ubi, lbj, ubj, n(ng), &
855 & imins, imaxs, jmins, jmaxs, &
856 & nnew, nstp, &
857 & gls, tke)
858!
859 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
860 CALL exchange_w3d_tile (ng, tile, &
861 & lbi, ubi, lbj, ubj, 0, n(ng), &
862 & tke(:,:,:,nnew))
863 CALL exchange_w3d_tile (ng, tile, &
864 & lbi, ubi, lbj, ubj, 0, n(ng), &
865 & gls(:,:,:,nnew))
866 CALL exchange_w3d_tile (ng, tile, &
867 & lbi, ubi, lbj, ubj, 0, n(ng), &
868 & akv)
869 DO itrc=1,nat
870 CALL exchange_w3d_tile (ng, tile, &
871 & lbi, ubi, lbj, ubj, 0, n(ng), &
872 & akt(:,:,:,itrc))
873 END DO
874 END IF
875
876# ifdef DISTRIBUTE
877 CALL mp_exchange3d (ng, tile, inlm, 3, &
878 & lbi, ubi, lbj, ubj, 0, n(ng), &
879 & nghostpoints, &
880 & ewperiodic(ng), nsperiodic(ng), &
881 & tke(:,:,:,nnew), &
882 & gls(:,:,:,nnew), &
883 & akv)
884 CALL mp_exchange4d (ng, tile, inlm, 1, &
885 & lbi, ubi, lbj, ubj, 0, n(ng), 1, nat, &
886 & nghostpoints, &
887 & ewperiodic(ng), nsperiodic(ng), &
888 & akt)
889# endif
890!
891 RETURN
892 END SUBROUTINE my25_corstep_tile
893#endif
894 END MODULE my25_corstep_mod
exchange_3d_mod
Definition exchange_3d.F:2
exchange_3d_mod::exchange_w3d_tile
subroutine exchange_w3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
Definition exchange_3d.F:920
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_mixing
Definition mod_mixing.F:2
mod_mixing::mixing
type(t_mixing), dimension(:), allocatable mixing
Definition mod_mixing.F:399
mod_ncparam
Definition mod_ncparam.F:2
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_param::lm
integer, dimension(:), allocatable lm
Definition mod_param.F:455
mod_param::mm
integer, dimension(:), allocatable mm
Definition mod_param.F:456
mod_scalars
Definition mod_scalars.F:2
mod_scalars::my_sh2
real(r8) my_sh2
Definition mod_scalars.F:1865
mod_scalars::my_sh1
real(r8) my_sh1
Definition mod_scalars.F:1864
mod_scalars::my_gh0
real(r8), parameter my_gh0
Definition mod_scalars.F:1855
mod_scalars::my_b1pm1o3
real(r8) my_b1pm1o3
Definition mod_scalars.F:1862
mod_scalars::my_sq
real(r8), parameter my_sq
Definition mod_scalars.F:1856
mod_scalars::my_lmax
real(r8), parameter my_lmax
Definition mod_scalars.F:1858
mod_scalars::vonkar
real(dp) vonkar
Definition mod_scalars.F:469
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::my_qmin
real(r8), parameter my_qmin
Definition mod_scalars.F:1859
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::iwest
integer, parameter iwest
Definition mod_scalars.F:1432
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::my_e1
real(r8), parameter my_e1
Definition mod_scalars.F:1853
mod_scalars::my_b1
real(r8), parameter my_b1
Definition mod_scalars.F:1848
mod_scalars::akk_bak
real(r8), dimension(:), allocatable akk_bak
Definition mod_scalars.F:889
mod_scalars::my_sm3
real(r8) my_sm3
Definition mod_scalars.F:1868
mod_scalars::my_sm4
real(r8) my_sm4
Definition mod_scalars.F:1869
mod_scalars::my_sm2
real(r8) my_sm2
Definition mod_scalars.F:1867
mod_scalars::akt_bak
real(r8), dimension(:,:), allocatable akt_bak
Definition mod_scalars.F:907
mod_scalars::akv_bak
real(r8), dimension(:), allocatable akv_bak
Definition mod_scalars.F:891
mod_scalars::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::my_b1p2o3
real(r8) my_b1p2o3
Definition mod_scalars.F:1861
mod_scalars::akt_limit
real(r8), dimension(:,:), allocatable akt_limit
Definition mod_scalars.F:908
mod_scalars::akv_limit
real(r8), dimension(:), allocatable akv_limit
Definition mod_scalars.F:892
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::my_e2
real(r8), parameter my_e2
Definition mod_scalars.F:1854
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange4d
subroutine mp_exchange4d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt, nghost, ew_periodic, ns_periodic, a, b, c)
Definition mp_exchange.F:3377
mp_exchange_mod::mp_exchange3d
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:2037
my25_corstep_mod
Definition my25_corstep.F:2
my25_corstep_mod::my25_corstep
subroutine, public my25_corstep(ng, tile)
Definition my25_corstep.F:36
my25_corstep_mod::my25_corstep_tile
subroutine my25_corstep_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, nnew, umask, vmask, huon, hvom, hz, pm, pn, z_r, z_w, u, v, w, bustr, bvstr, sustr, svstr, bvf, akt, akv, akk, lscale, gls, tke)
Definition my25_corstep.F:108
tkebc_mod
Definition tkebc_im.F:2
tkebc_mod::tkebc_tile
subroutine, public tkebc_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, nout, nstp, gls, tke)
Definition tkebc_im.F:56
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3