ROMS
Loading...
Searching...
No Matches
gls_corstep.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE gls_corstep_mod
3#if defined NONLINEAR && defined GLS_MIXING && defined SOLVE3D
4!
5!git $Id$
6!=======================================================================
7! Copyright (c) 2002-2025 The ROMS Group !
8! Licensed under a MIT/X style license Hernan G. Arango !
9! See License_ROMS.md Alexander F. Shchepetkin !
10!==================================================== John C. Warner ===
11! !
12! This routine perfoms the corrector step for turbulent kinetic !
13! energy and generic length scale prognostic variables, tke and !
14! gls. !
15! !
16! References: !
17! !
18! Umlauf, L. and H. Burchard, 2003: A generic length-scale !
19! Equation for geophysical turbulence models, J. Marine Res., !
20! 61, 235-265. !
21! !
22! Warner, J. C., C. R. Sherwood, H. G. Arango and R. P. Signell, !
23! 2005: Performance of four turbulence closure models !
24! implemented using a generic length scale method, Ocean !
25! Modelling, 8, 81-113. !
26! !
27!=======================================================================
28!
29 implicit none
30!
31 PRIVATE
32 PUBLIC :: gls_corstep
33!
34 CONTAINS
35!
36!***********************************************************************
37 SUBROUTINE gls_corstep (ng, tile)
38!***********************************************************************
39!
40 USE mod_param
41 USE mod_forces
42 USE mod_grid
43 USE mod_mixing
44 USE mod_ocean
45 USE mod_stepping
46!
47! Imported variable declarations.
48!
49 integer, intent(in) :: ng, tile
50!
51! Local variable declarations.
52!
53 character (len=*), parameter :: myfile = &
54 & __FILE__
55!
56# include "tile.h"
57!
58# ifdef PROFILE
59 CALL wclock_on (ng, inlm, 19, __line__, myfile)
60# endif
61 CALL gls_corstep_tile (ng, tile, &
62 & lbi, ubi, lbj, ubj, &
63 & imins, imaxs, jmins, jmaxs, &
64 & nstp(ng), nnew(ng), &
65# ifdef MASKING
66 & grid(ng) % umask, &
67 & grid(ng) % vmask, &
68# endif
69 & grid(ng) % Huon, &
70 & grid(ng) % Hvom, &
71 & grid(ng) % Hz, &
72 & grid(ng) % pm, &
73 & grid(ng) % pn, &
74 & grid(ng) % z_r, &
75 & grid(ng) % z_w, &
76 & grid(ng) % ZoBot, &
77 & ocean(ng) % u, &
78 & ocean(ng) % v, &
79 & ocean(ng) % W, &
80# ifdef WEC_VF
81 & ocean(ng) % W_stokes, &
82# endif
83 & forces(ng) % bustr, &
84 & forces(ng) % bvstr, &
85 & forces(ng) % sustr, &
86 & forces(ng) % svstr, &
87# ifdef ZOS_HSIG
88 & forces(ng) % Hwave, &
89# endif
90# ifdef TKE_WAVEDISS
91 & forces(ng) % Dissip_break, &
92 & forces(ng) % Dissip_wcap, &
93# endif
94 & mixing(ng) % Akt, &
95 & mixing(ng) % Akv, &
96 & mixing(ng) % bvf, &
97 & mixing(ng) % Akk, &
98 & mixing(ng) % Akp, &
99 & mixing(ng) % Lscale, &
100 & mixing(ng) % gls, &
101 & mixing(ng) % tke)
102# ifdef PROFILE
103 CALL wclock_off (ng, inlm, 19, __line__, myfile)
104# endif
105!
106 RETURN
107 END SUBROUTINE gls_corstep
108!
109!***********************************************************************
110 SUBROUTINE gls_corstep_tile (ng, tile, &
111 & LBi, UBi, LBj, UBj, &
112 & IminS, ImaxS, JminS, JmaxS, &
113 & nstp, nnew, &
114# ifdef MASKING
115 & umask, vmask, &
116# endif
117 & Huon, Hvom, Hz, pm, pn, z_r, z_w, &
118 & ZoBot, &
119 & u, v, W, &
120# ifdef WEC_VF
121 & W_stokes, &
122# endif
123 & bustr, bvstr, sustr, svstr, &
124# ifdef ZOS_HSIG
125 & Hwave, &
126# endif
127# ifdef TKE_WAVEDISS
128 & Dissip_break, Dissip_wcap, &
129# endif
130 & Akt, Akv, bvf, &
131 & Akk, Akp, Lscale, gls, tke)
132!***********************************************************************
133!
134 USE mod_param
135 USE mod_ncparam
136 USE mod_scalars
137!
138 USE exchange_3d_mod, ONLY : exchange_w3d_tile
139# ifdef DISTRIBUTE
140 USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
141# endif
142 USE tkebc_mod, ONLY : tkebc_tile
143!
144! Imported variable declarations.
145!
146 integer, intent(in) :: ng, tile
147 integer, intent(in) :: LBi, UBi, LBj, UBj
148 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
149 integer, intent(in) :: nstp, nnew
150!
151# ifdef ASSUMED_SHAPE
152# ifdef MASKING
153 real(r8), intent(in) :: umask(LBi:,LBj:)
154 real(r8), intent(in) :: vmask(LBi:,LBj:)
155# endif
156 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
157 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
158 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
159 real(r8), intent(in) :: pm(LBi:,LBj:)
160 real(r8), intent(in) :: pn(LBi:,LBj:)
161 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
162 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
163 real(r8), intent(in) :: ZoBot(LBi:,LBj:)
164 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
165 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
166 real(r8), intent(in) :: W(LBi:,LBj:,0:)
167# ifdef WEC_VF
168 real(r8), intent(in) :: W_stokes(LBi:,LBj:,0:)
169# endif
170 real(r8), intent(in) :: bustr(LBi:,LBj:)
171 real(r8), intent(in) :: bvstr(LBi:,LBj:)
172 real(r8), intent(in) :: sustr(LBi:,LBj:)
173 real(r8), intent(in) :: svstr(LBi:,LBj:)
174# ifdef ZOS_HSIG
175 real(r8), intent(in) :: Hwave(LBi:,LBj:)
176# endif
177# ifdef TKE_WAVEDISS
178 real(r8), intent(in) :: Dissip_break(LBi:,LBj:)
179 real(r8), intent(in) :: Dissip_wcap(LBi:,LBj:)
180# endif
181 real(r8), intent(in) :: bvf(LBi:,LBj:,0:)
182
183 real(r8), intent(inout) :: Akt(LBi:,LBj:,0:,:)
184 real(r8), intent(inout) :: Akv(LBi:,LBj:,0:)
185 real(r8), intent(inout) :: Akk(LBi:,LBj:,0:)
186 real(r8), intent(inout) :: Akp(LBi:,LBj:,0:)
187 real(r8), intent(inout) :: Lscale(LBi:,LBj:,0:)
188 real(r8), intent(inout) :: gls(LBi:,LBj:,0:,:)
189 real(r8), intent(inout) :: tke(LBi:,LBj:,0:,:)
190# else
191# ifdef MASKING
192 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
193 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
194# endif
195 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
196 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
197 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
198 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
199 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
200 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
201 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
202 real(r8), intent(in) :: ZoBot(LBi:UBi,LBj:UBj)
203 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
204 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
205 real(r8), intent(in) :: W(LBi:UBi,LBj:UBj,0:N(ng))
206# ifdef WEC_VF
207 real(r8), intent(in) :: W_stokes(LBi:UBi,LBj:UBj,0:N(ng))
208# endif
209 real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj)
210 real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
211 real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj)
212 real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj)
213# if defined ZOS_HSIG
214 real(r8), intent(in) :: Hwave(LBi:UBi,LBj:UBj)
215# endif
216# ifdef TKE_WAVEDISS
217 real(r8), intent(in) :: Dissip_break(LBi:UBi,LBj:UBj)
218 real(r8), intent(in) :: Dissip_wcap(LBi:UBi,LBj:UBj)
219# endif
220 real(r8), intent(in) :: bvf(LBi:UBi,LBj:UBj,0:N(ng))
221
222 real(r8), intent(inout) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
223 real(r8), intent(inout) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
224 real(r8), intent(inout) :: Akk(LBi:UBi,LBj:UBj,0:N(ng))
225 real(r8), intent(inout) :: Akp(LBi:UBi,LBj:UBj,0:N(ng))
226 real(r8), intent(inout) :: Lscale(LBi:UBi,LBj:UBj,0:N(ng))
227 real(r8), intent(inout) :: gls(LBi:UBi,LBj:UBj,0:N(ng),3)
228 real(r8), intent(inout) :: tke(LBi:UBi,LBj:UBj,0:N(ng),3)
229# endif
230!
231! Local variable declarations.
232!
233 logical :: Lmy25
234 integer :: i, itrc, j, k
235
236 real(r8), parameter :: Gadv = 1.0_r8/3.0_r8
237 real(r8), parameter :: eps = 1.0e-10_r8
238 real(r8) :: Zos_min
239
240 real(r8) :: Gh, Gm, Kprod, Ls_unlmt, Ls_lmt, Pprod, Sh, Sm
241 real(r8) :: cff, cff1, cff2, cff3
242 real(r8) :: cmu_fac1, cmu_fac2, cmu_fac3, cmu_fac4
243 real(r8) :: gls_c3, gls_exp1, gls_fac1, gls_fac2, gls_fac3
244 real(r8) :: gls_fac4, gls_fac5, gls_fac6, ql, sqrt2, strat2
245 real(r8) :: tke_exp1, tke_exp2, tke_exp3, tke_exp4, wall_fac
246 real(r8) :: gls_d, gls_sigp_cb, ogls_sigp, sig_eff
247 real(r8) :: L_sft
248# if defined CRAIG_BANNER || defined TKE_WAVEDISS
249 real(r8) :: cb_wallE
250# endif
251
252 real(r8), dimension(IminS:ImaxS) :: tke_fluxt
253 real(r8), dimension(IminS:ImaxS) :: tke_fluxb
254 real(r8), dimension(IminS:ImaxS) :: gls_fluxt
255 real(r8), dimension(IminS:ImaxS) :: gls_fluxb
256 real(r8), dimension(IminS:ImaxS) :: Zos_eff
257
258 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BCK
259 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BCP
260 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
261 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCK
262 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCP
263 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dU
264 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dV
265
266 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: shear2
267 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: buoy2
268
269 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FEK
270 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FEP
271 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FXK
272 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FXP
273 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Zob_min
274 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curvK
275 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curvP
276 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradK
277 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradP
278
279# include "set_bounds.h"
280!
281!-----------------------------------------------------------------------
282! Compute several constants.
283!-----------------------------------------------------------------------
284!
285 zos_min=max(zos(ng),0.0001_r8)
286 DO j=jstr,jend
287 DO i=istr,iend
288 zob_min(i,j)=max(zobot(i,j),0.0001_r8)
289 END DO
290 END DO
291!
292 lmy25=.false.
293 IF ((gls_p(ng).eq.0.0_r8).and. &
294 & (gls_n(ng).eq.1.0_r8).and. &
295 & (gls_m(ng).eq.1.0_r8)) THEN
296 lmy25=.true.
297 END IF
298
299# if defined CRAIG_BANNER || defined TKE_WAVEDISS
300 IF (lmy25) THEN
301!! cb_wallE=1.0_r8+gls_E2
302 cb_walle=1.25_r8
303 ELSE
304 cb_walle=1.0_r8
305 END IF
306 l_sft=vonkar
307 cff1=sqrt(1.5_r8*gls_sigk(ng))*gls_cmu0(ng)/l_sft
308 gls_sigp_cb=l_sft**2/(gls_cmu0(ng)**2*gls_c2(ng)*cb_walle)* &
309 & (gls_n(ng)**2- &
310 & cff1*gls_n(ng)/3.0_r8*(4.0_r8*gls_m(ng)+1.0_r8)+ &
311 & cff1**2*gls_m(ng)/9.0_r8*(2.0_r8+4.0_r8*gls_m(ng)))
312# else
313 l_sft=vonkar
314 gls_sigp_cb=gls_sigp(ng)
315# endif
316 ogls_sigp=1.0_r8/gls_sigp_cb
317!
318 sqrt2=sqrt(2.0_r8)
319 cmu_fac1=gls_cmu0(ng)**(-gls_p(ng)/gls_n(ng))
320 cmu_fac2=gls_cmu0(ng)**(3.0_r8+gls_p(ng)/gls_n(ng))
321 cmu_fac3=1.0_r8/gls_cmu0(ng)**2.0_r8
322 cmu_fac4=(1.5_r8*gls_sigk(ng))**(1.0_r8/3.0_r8)/ &
323 & (gls_cmu0(ng)**(4.0_r8/3.0_r8))
324 gls_fac1=gls_n(ng)*gls_cmu0(ng)**(gls_p(ng)+1.0_r8)
325 gls_fac2=gls_cmu0(ng)**(gls_p(ng))*gls_n(ng)* &
326 & vonkar**(gls_n(ng))
327 gls_fac3=gls_cmu0(ng)**(gls_p(ng))*gls_n(ng)
328 gls_fac4=gls_cmu0(ng)**(gls_p(ng))
329 gls_fac5=0.56_r8**(0.5_r8*gls_n(ng))*gls_cmu0(ng)**gls_p(ng)
330 gls_fac6=8.0_r8/gls_cmu0(ng)**6.0_r8
331!
332 gls_exp1=1.0_r8/gls_n(ng)
333 tke_exp1=gls_m(ng)/gls_n(ng)
334 tke_exp2=0.5_r8+gls_m(ng)/gls_n(ng)
335 tke_exp3=0.5_r8+gls_m(ng)
336 tke_exp4=gls_m(ng)+0.5_r8*gls_n(ng)
337!
338!-----------------------------------------------------------------------
339! Compute vertical velocity shear at W-points.
340!-----------------------------------------------------------------------
341!
342# ifdef RI_SPLINES
343 DO j=jstrm1,jendp1
344 DO i=istrm1,iendp1
345 cf(i,0)=0.0_r8
346 du(i,0)=0.0_r8
347 dv(i,0)=0.0_r8
348 END DO
349 DO k=1,n(ng)-1
350 DO i=istrm1,iendp1
351 cff=1.0_r8/(2.0_r8*hz(i,j,k+1)+ &
352 & hz(i,j,k)*(2.0_r8-cf(i,k-1)))
353 cf(i,k)=cff*hz(i,j,k+1)
354 du(i,k)=cff*(3.0_r8*(u(i ,j,k+1,nstp)-u(i, j,k,nstp)+ &
355 & u(i+1,j,k+1,nstp)-u(i+1,j,k,nstp))- &
356 & hz(i,j,k)*du(i,k-1))
357 dv(i,k)=cff*(3.0_r8*(v(i,j ,k+1,nstp)-v(i,j ,k,nstp)+ &
358 & v(i,j+1,k+1,nstp)-v(i,j+1,k,nstp))- &
359 & hz(i,j,k)*dv(i,k-1))
360 END DO
361 END DO
362 DO i=istrm1,iendp1
363 du(i,n(ng))=0.0_r8
364 dv(i,n(ng))=0.0_r8
365 END DO
366 DO k=n(ng)-1,1,-1
367 DO i=istrm1,iendp1
368 du(i,k)=du(i,k)-cf(i,k)*du(i,k+1)
369 dv(i,k)=dv(i,k)-cf(i,k)*dv(i,k+1)
370 END DO
371 END DO
372 DO k=1,n(ng)-1
373 DO i=istrm1,iendp1
374 shear2(i,j,k)=du(i,k)*du(i,k)+dv(i,k)*dv(i,k)
375 END DO
376 END DO
377 END DO
378# else
379 DO k=1,n(ng)-1
380 DO j=jstrm1,jendp1
381 DO i=istrm1,iendp1
382 cff=0.5_r8/(z_r(i,j,k+1)-z_r(i,j,k))
383 shear2(i,j,k)=(cff*(u(i ,j,k+1,nstp)-u(i ,j,k,nstp)+ &
384 & u(i+1,j,k+1,nstp)-u(i+1,j,k,nstp)))**2+ &
385 & (cff*(v(i,j ,k+1,nstp)-v(i,j ,k,nstp)+ &
386 & v(i,j+1,k+1,nstp)-v(i,j+1,k,nstp)))**2
387 END DO
388 END DO
389 END DO
390# endif
391!
392! Load Brunt-Vaisala frequency.
393!
394 DO k=1,n(ng)-1
395 DO j=jstr-1,jend+1
396 DO i=istr-1,iend+1
397 buoy2(i,j,k)=bvf(i,j,k)
398 END DO
399 END DO
400 END DO
401# ifdef N2S2_HORAVG
402!
403!-----------------------------------------------------------------------
404! Smooth horizontally buoyancy and shear. Use buoy2(:,:,0) and
405! shear2(:,:,0) as scratch utility array.
406!-----------------------------------------------------------------------
407!
408 DO k=1,n(ng)-1
409 IF (domain(ng)%Western_Edge(tile)) THEN
410 DO j=max(1,jstr-1),min(jend+1,mm(ng))
411 shear2(istr-1,j,k)=shear2(istr,j,k)
412 END DO
413 END IF
414 IF (domain(ng)%Eastern_Edge(tile)) THEN
415 DO j=max(1,jstr-1),min(jend+1,mm(ng))
416 shear2(iend+1,j,k)=shear2(iend,j,k)
417 END DO
418 END IF
419 IF (domain(ng)%Southern_Edge(tile)) THEN
420 DO i=max(1,istr-1),min(iend+1,lm(ng))
421 shear2(i,jstr-1,k)=shear2(i,jstr,k)
422 END DO
423 END IF
424 IF (domain(ng)%Northern_Edge(tile)) THEN
425 DO i=max(1,istr-1),min(iend+1,lm(ng))
426 shear2(i,jend+1,k)=shear2(i,jend,k)
427 END DO
428 END IF
429 IF (domain(ng)%SouthWest_Corner(tile)) THEN
430 shear2(istr-1,jstr-1,k)=shear2(istr,jstr,k)
431 END IF
432 IF (domain(ng)%NorthWest_Corner(tile)) THEN
433 shear2(istr-1,jend+1,k)=shear2(istr,jend,k)
434 END IF
435 IF (domain(ng)%SouthEast_Corner(tile)) THEN
436 shear2(iend+1,jstr-1,k)=shear2(iend,jstr,k)
437 END IF
438 IF (domain(ng)%NorthEast_Corner(tile)) THEN
439 shear2(iend+1,jend+1,k)=shear2(iend,jend,k)
440 END IF
441!
442! Average horizontally.
443!
444 DO j=jstr-1,jend
445 DO i=istr-1,iend
446 buoy2(i,j,0)=0.25_r8*(buoy2(i,j ,k)+buoy2(i+1,j ,k)+ &
447 & buoy2(i,j+1,k)+buoy2(i+1,j+1,k))
448 shear2(i,j,0)=0.25_r8*(shear2(i,j ,k)+shear2(i+1,j ,k)+ &
449 & shear2(i,j+1,k)+shear2(i+1,j+1,k))
450 END DO
451 END DO
452 DO j=jstr,jend
453 DO i=istr,iend
454 buoy2(i,j,k)=0.25_r8*(buoy2(i,j ,0)+buoy2(i-1,j ,0)+ &
455 & buoy2(i,j-1,0)+buoy2(i-1,j-1,0))
456 shear2(i,j,k)=0.25_r8*(shear2(i,j ,0)+shear2(i-1,j ,0)+ &
457 & shear2(i,j-1,0)+shear2(i-1,j-1,0))
458 END DO
459 END DO
460 END DO
461# endif
462!
463!-----------------------------------------------------------------------
464! Time-step advective terms.
465!-----------------------------------------------------------------------
466!
467! At entry, it is assumed that the turbulent kinetic energy fields
468! "tke" and "gls", at time level "nnew", are set to its values at
469! time level "nstp" multiplied by the grid box thicknesses Hz
470! (from old time step and at W-points).
471!
472 DO k=1,n(ng)-1
473# ifdef K_C2ADVECTION
474!
475! Second-order, centered differences advection.
476!
477 DO j=jstr,jend
478 DO i=istr,iend+1
479 cff=0.25_r8*(huon(i,j,k)+huon(i,j,k+1))
480 fxk(i,j)=cff*(tke(i,j,k,3)+tke(i-1,j,k,3))
481 fxp(i,j)=cff*(gls(i,j,k,3)+gls(i-1,j,k,3))
482 END DO
483 END DO
484 DO j=jstr,jend+1
485 DO i=istr,iend
486 cff=0.25_r8*(hvom(i,j,k)+hvom(i,j,k+1))
487 fek(i,j)=cff*(tke(i,j,k,3)+tke(i,j-1,k,3))
488 fep(i,j)=cff*(gls(i,j,k,3)+gls(i,j-1,k,3))
489 END DO
490 END DO
491# else
492 DO j=jstr,jend
493 DO i=istrm1,iendp2
494 gradk(i,j)=(tke(i,j,k,3)-tke(i-1,j,k,3))
495# ifdef MASKING
496 gradk(i,j)=gradk(i,j)*umask(i,j)
497# endif
498 gradp(i,j)=(gls(i,j,k,3)-gls(i-1,j,k,3))
499# ifdef MASKING
500 gradp(i,j)=gradp(i,j)*umask(i,j)
501# endif
502 END DO
503 END DO
504 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
505 IF (domain(ng)%Western_Edge(tile)) THEN
506 DO j=jstr,jend
507 gradk(istr-1,j)=gradk(istr,j)
508 gradp(istr-1,j)=gradp(istr,j)
509 END DO
510 END IF
511 END IF
512 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
513 IF (domain(ng)%Eastern_Edge(tile)) THEN
514 DO j=jstr,jend
515 gradk(iend+2,j)=gradk(iend+1,j)
516 gradp(iend+2,j)=gradp(iend+1,j)
517 END DO
518 END IF
519 END IF
520# ifdef K_C4ADVECTION
521!
522! Fourth-order, centered differences advection.
523!
524 cff1=1.0_r8/6.0_r8
525 DO j=jstr,jend
526 DO i=istr,iend+1
527 cff=0.5_r8*(huon(i,j,k)+huon(i,j,k+1))
528 fxk(i,j)=cff*0.5_r8*(tke(i-1,j,k,3)+tke(i,j,k,3)- &
529 & cff1*(gradk(i+1,j)-gradk(i-1,j)))
530 fxp(i,j)=cff*0.5_r8*(gls(i-1,j,k,3)+gls(i,j,k,3)- &
531 & cff1*(gradp(i+1,j)-gradp(i-1,j)))
532 END DO
533 END DO
534# else
535!
536! Third-order, upstream bias advection with velocity dependent
537! hyperdiffusion.
538!
539 DO j=jstr,jend
540 DO i=istr-1,iend+1
541 curvk(i,j)=gradk(i+1,j)-gradk(i,j)
542 curvp(i,j)=gradp(i+1,j)-gradp(i,j)
543 END DO
544 END DO
545 DO j=jstr,jend
546 DO i=istr,iend+1
547 cff=0.5_r8*(huon(i,j,k)+huon(i,j,k+1))
548 IF (cff.gt.0.0_r8) THEN
549 cff1=curvk(i-1,j)
550 cff2=curvp(i-1,j)
551 ELSE
552 cff1=curvk(i,j)
553 cff2=curvp(i,j)
554 END IF
555 fxk(i,j)=cff*0.5_r8*(tke(i-1,j,k,3)+tke(i,j,k,3)- &
556 & gadv*cff1)
557 fxp(i,j)=cff*0.5_r8*(gls(i-1,j,k,3)+gls(i,j,k,3)- &
558 & gadv*cff2)
559 END DO
560 END DO
561# endif
562 DO j=jstrm1,jendp2
563 DO i=istr,iend
564 gradk(i,j)=(tke(i,j,k,3)-tke(i,j-1,k,3))
565# ifdef MASKING
566 gradk(i,j)=gradk(i,j)*vmask(i,j)
567# endif
568 gradp(i,j)=(gls(i,j,k,3)-gls(i,j-1,k,3))
569# ifdef MASKING
570 gradp(i,j)=gradp(i,j)*vmask(i,j)
571# endif
572 END DO
573 END DO
574 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
575 IF (domain(ng)%Southern_Edge(tile)) THEN
576 DO i=istr,iend
577 gradk(i,jstr-1)=gradk(i,jstr)
578 gradp(i,jstr-1)=gradp(i,jstr)
579 END DO
580 END IF
581 END IF
582 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
583 IF (domain(ng)%Northern_Edge(tile)) THEN
584 DO i=istr,iend
585 gradk(i,jend+2)=gradk(i,jend+1)
586 gradp(i,jend+2)=gradp(i,jend+1)
587 END DO
588 END IF
589 END IF
590# ifdef K_C4ADVECTION
591 cff1=1.0_r8/6.0_r8
592 DO j=jstr,jend+1
593 DO i=istr,iend
594 cff=0.5_r8*(hvom(i,j,k)+hvom(i,j,k+1))
595 fek(i,j)=cff*0.5_r8*(tke(i,j-1,k,3)+tke(i,j,k,3)- &
596 & cff1*(gradk(i,j+1)-gradk(i,j-1)))
597 fep(i,j)=cff*0.5_r8*(gls(i,j-1,k,3)+gls(i,j,k,3)- &
598 & cff1*(gradp(i,j+1)-gradp(i,j-1)))
599 END DO
600 END DO
601# else
602 DO j=jstr-1,jend+1
603 DO i=istr,iend
604 curvk(i,j)=gradk(i,j+1)-gradk(i,j)
605 curvp(i,j)=gradp(i,j+1)-gradp(i,j)
606 END DO
607 END DO
608 DO j=jstr,jend+1
609 DO i=istr,iend
610 cff=0.5_r8*(hvom(i,j,k)+hvom(i,j,k+1))
611 IF (cff.gt.0.0_r8) THEN
612 cff1=curvk(i,j-1)
613 cff2=curvp(i,j-1)
614 ELSE
615 cff1=curvk(i,j)
616 cff2=curvp(i,j)
617 END IF
618 fek(i,j)=cff*0.5_r8*(tke(i,j-1,k,3)+tke(i,j,k,3)- &
619 & gadv*cff1)
620 fep(i,j)=cff*0.5_r8*(gls(i,j-1,k,3)+gls(i,j,k,3)- &
621 & gadv*cff2)
622 END DO
623 END DO
624# endif
625# endif
626!
627! Time-step horizontal advection.
628!
629 DO j=jstr,jend
630 DO i=istr,iend
631 cff=dt(ng)*pm(i,j)*pn(i,j)
632 tke(i,j,k,nnew)=tke(i,j,k,nnew)- &
633 & cff*(fxk(i+1,j)-fxk(i,j)+ &
634 & fek(i,j+1)-fek(i,j))
635 tke(i,j,k,nnew)=max(tke(i,j,k,nnew),gls_kmin(ng))
636 gls(i,j,k,nnew)=gls(i,j,k,nnew)- &
637 & cff*(fxp(i+1,j)-fxp(i,j)+ &
638 & fep(i,j+1)-fep(i,j))
639 gls(i,j,k,nnew)=max(gls(i,j,k,nnew),gls_pmin(ng))
640 END DO
641 END DO
642 END DO
643!
644! Compute vertical advection.
645!
646 DO j=jstr,jend
647# ifdef K_C2ADVECTION
648 DO k=1,n(ng)
649 DO i=istr,iend
650 cff=0.25_r8*(w(i,j,k)+w(i,j,k-1))
651# ifdef WEC_VF
652 cff=cff+0.25_r8*(w_stokes(i,j,k)+w_stokes(i,j,k-1))
653# endif
654 fck(i,k)=cff*(tke(i,j,k,3)+tke(i,j,k-1,3))
655 fcp(i,k)=cff*(gls(i,j,k,3)+gls(i,j,k-1,3))
656 END DO
657 END DO
658# else
659 cff1=7.0_r8/12.0_r8
660 cff2=1.0_r8/12.0_r8
661 DO k=2,n(ng)-1
662 DO i=istr,iend
663 cff=0.5*(w(i,j,k)+w(i,j,k-1))
664# ifdef WEC_VF
665 cff=cff+0.5_r8*(w_stokes(i,j,k)+w_stokes(i,j,k-1))
666# endif
667 fck(i,k)=cff*(cff1*(tke(i,j,k-1,3)+ &
668 & tke(i,j,k ,3))- &
669 & cff2*(tke(i,j,k-2,3)+ &
670 & tke(i,j,k+1,3)))
671 fcp(i,k)=cff*(cff1*(gls(i,j,k-1,3)+ &
672 & gls(i,j,k ,3))- &
673 & cff2*(gls(i,j,k-2,3)+ &
674 & gls(i,j,k+1,3)))
675 END DO
676 END DO
677 cff1=1.0_r8/3.0_r8
678 cff2=5.0_r8/6.0_r8
679 cff3=1.0_r8/6.0_r8
680 DO i=istr,iend
681 cff=0.5_r8*(w(i,j,0)+w(i,j,1))
682# ifdef WEC_VF
683 cff=cff+0.5_r8*(w_stokes(i,j,0)+w_stokes(i,j,1))
684# endif
685 fck(i,1)=cff*(cff1*tke(i,j,0,3)+ &
686 & cff2*tke(i,j,1,3)- &
687 & cff3*tke(i,j,2,3))
688 fcp(i,1)=cff*(cff1*gls(i,j,0,3)+ &
689 & cff2*gls(i,j,1,3)- &
690 & cff3*gls(i,j,2,3))
691 cff=0.5_r8*(w(i,j,n(ng))+w(i,j,n(ng)-1))
692# ifdef WEC_VF
693 cff=cff+0.5_r8*(w_stokes(i,j,n(ng))+w_stokes(i,j,n(ng)-1))
694# endif
695 fck(i,n(ng))=cff*(cff1*tke(i,j,n(ng) ,3)+ &
696 & cff2*tke(i,j,n(ng)-1,3)- &
697 & cff3*tke(i,j,n(ng)-2,3))
698 fcp(i,n(ng))=cff*(cff1*gls(i,j,n(ng) ,3)+ &
699 & cff2*gls(i,j,n(ng)-1,3)- &
700 & cff3*gls(i,j,n(ng)-2,3))
701 END DO
702# endif
703!
704! Time-step vertical advection term.
705!
706 DO k=1,n(ng)-1
707 DO i=istr,iend
708 cff=dt(ng)*pm(i,j)*pn(i,j)
709 tke(i,j,k,nnew)=tke(i,j,k,nnew)- &
710 & cff*(fck(i,k+1)-fck(i,k))
711 tke(i,j,k,nnew)=max(tke(i,j,k,nnew),gls_kmin(ng))
712 gls(i,j,k,nnew)=gls(i,j,k,nnew)- &
713 & cff*(fcp(i,k+1)-fcp(i,k))
714 gls(i,j,k,nnew)=max(gls(i,j,k,nnew),gls_pmin(ng))
715 END DO
716 END DO
717!
718!----------------------------------------------------------------------
719! Compute vertical mixing, turbulent production and turbulent
720! dissipation terms.
721!----------------------------------------------------------------------
722!
723! Set term for vertical mixing of turbulent fields.
724!
725 cff=-0.5_r8*dt(ng)
726 DO i=istr,iend
727 DO k=2,n(ng)-1
728 fck(i,k)=cff*(akk(i,j,k)+akk(i,j,k-1))/hz(i,j,k)
729 fcp(i,k)=cff*(akp(i,j,k)+akp(i,j,k-1))/hz(i,j,k)
730 cf(i,k)=0.0_r8
731 END DO
732 fcp(i,1)=0.0_r8
733 fcp(i,n(ng))=0.0_r8
734 fck(i,1)=0.0_r8
735 fck(i,n(ng))=0.0_r8
736 END DO
737!
738! Compute production and dissipation terms.
739!
740 DO i=istr,iend
741 DO k=1,n(ng)-1
742!
743! Compute shear and bouyant production of turbulent energy (m3/s3)
744! at W-points (ignore small negative values of buoyancy).
745!
746 strat2=buoy2(i,j,k)
747 IF (strat2.gt.0.0_r8) THEN
748 gls_c3=gls_c3m(ng)
749 ELSE
750 gls_c3=gls_c3p(ng)
751 END IF
752 kprod=shear2(i,j,k)*(akv(i,j,k)-akv_bak(ng))- &
753 & strat2*(akt(i,j,k,itemp)-akt_bak(itemp,ng))
754 pprod=gls_c1(ng)*shear2(i,j,k)*(akv(i,j,k)-akv_bak(ng))- &
755 & gls_c3*strat2*(akt(i,j,k,itemp)-akt_bak(itemp,ng))
756!
757! If negative production terms, then add buoyancy to dissipation terms
758! (BCK and BCP) below, using "cff1" and "cff2" as the on/off switch.
759!
760 cff1=1.0_r8
761 IF (kprod.lt.0.0_r8) THEN
762 kprod=kprod+strat2*(akt(i,j,k,itemp)-akt_bak(itemp,ng))
763 cff1=0.0_r8
764 END IF
765 cff2=1.0_r8
766 IF (pprod.lt.0.0_r8) THEN
767 pprod=pprod+gls_c3*strat2*(akt(i,j,k,itemp)- &
768 & akt_bak(itemp,ng))
769 cff2=0.0_r8
770 END IF
771!
772! Time-step shear and buoyancy production terms.
773!
774 cff=0.5_r8*(hz(i,j,k)+hz(i,j,k+1))
775 tke(i,j,k,nnew)=tke(i,j,k,nnew)+ &
776 & dt(ng)*cff*kprod
777 gls(i,j,k,nnew)=gls(i,j,k,nnew)+ &
778 & dt(ng)*cff*pprod*gls(i,j,k,nstp)/ &
779 & max(tke(i,j,k,nstp),gls_kmin(ng))
780!
781! Compute dissipation of turbulent energy (m3/s3).
782!
783 wall_fac=1.0_r8
784 IF (lmy25) THEN
785!
786! Parabolic wall function, L = ds db / (ds + db).
787!
788!! wall_fac=1.0_r8+gls_E2/(vonKar*vonKar)* &
789!! & (gls(i,j,k,nstp)**( gls_exp1)*cmu_fac1* &
790!! & tke(i,j,k,nstp)**(-tke_exp1)* &
791!! & (1.0_r8/(z_w(i,j,N(ng))-z_w(i,j,k))+ &
792!! & 1.0_r8/(z_w(i,j,k)-z_w(i,j,0))))**2
793!!
794!! Triangular wall function, L = min (ds, db).
795!!
796!! wall_fac=1.0_r8+gls_E2/(vonKar*vonKar)* &
797!! & (gls(i,j,k,nstp)**( gls_exp1)*cmu_fac1* &
798!! & tke(i,j,k,nstp)**(-tke_exp1)* &
799!! & (1.0_r8/MIN((z_w(i,j,N(ng))-z_w(i,j,k)), &
800!! & (z_w(i,j,k)-z_w(i,j,0)))))**2
801!!
802!! Linear wall function for , L = ds (=dist to surface).
803!!
804!! wall_fac=1.0_r8+gls_E2/(vonKar*vonKar)* &
805!! & (gls(i,j,k,nstp)**( gls_exp1)*cmu_fac1* &
806!! & tke(i,j,k,nstp)**(-tke_exp1)* &
807!! & (1.0_r8/ (z_w(i,j,N(ng))-z_w(i,j,k))))**2
808!!
809!! Parabolic wall function + free surface correction
810!
811 wall_fac=1.0_r8+gls_e2/(vonkar*vonkar)* &
812 & (gls(i,j,k,nstp)**( gls_exp1)*cmu_fac1* &
813 & tke(i,j,k,nstp)**(-tke_exp1)* &
814 & (1.0_r8/ (z_w(i,j,k)-z_w(i,j,0))))**2+ &
815 & 0.25_r8/(vonkar*vonkar)* &
816 & (gls(i,j,k,nstp)**( gls_exp1)*cmu_fac1* &
817 & tke(i,j,k,nstp)**(-tke_exp1)* &
818 & (1.0_r8/ (z_w(i,j,n(ng))-z_w(i,j,k))))**2
819 END IF
820!
821 bck(i,k)=cff*(1.0_r8+dt(ng)* &
822 & gls(i,j,k,nstp)**(-gls_exp1)*cmu_fac2* &
823 & tke(i,j,k,nstp)**( tke_exp2)+ &
824 & dt(ng)*(1.0_r8-cff1)*strat2* &
825 & (akt(i,j,k,itemp)-akt_bak(itemp,ng))/ &
826 & tke(i,j,k,nstp))- &
827 & fck(i,k)-fck(i,k+1)
828 bcp(i,k)=cff*(1.0_r8+dt(ng)*gls_c2(ng)*wall_fac* &
829 & gls(i,j,k,nstp)**(-gls_exp1)*cmu_fac2* &
830 & tke(i,j,k,nstp)**( tke_exp2)+ &
831 & dt(ng)*(1.0_r8-cff2)*gls_c3*strat2* &
832 & (akt(i,j,k,itemp)-akt_bak(itemp,ng))/ &
833 & tke(i,j,k,nstp))- &
834 & fcp(i,k)-fcp(i,k+1)
835 END DO
836 END DO
837!
838!-----------------------------------------------------------------------
839! Time-step dissipation and vertical diffusion terms implicitly.
840!-----------------------------------------------------------------------
841!
842! Set Dirichlet surface and bottom boundary conditions. Compute
843! surface roughness from wind stress (Charnock) and set Craig and
844! Banner wave breaking surface flux, if appropriate.
845!
846 DO i=istr,iend
847# if defined CRAIG_BANNER
848 tke(i,j,n(ng),nnew)=max(cmu_fac4*0.5_r8* &
849 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
850 & (svstr(i,j)+svstr(i,j+1))**2)* &
851 & crgban_cw(ng)**(2.0_r8/3.0_r8), &
852 & gls_kmin(ng))
853# elif defined TKE_WAVEDISS
854 tke(i,j,n(ng),nnew)=max(cmu_fac4*(sz_alpha(ng)* &
855 & (dissip_break(i,j)+ &
856 & dissip_wcap(i,j)))** &
857 & (2.0_r8/3.0_r8), gls_kmin(ng))
858# else
859 tke(i,j,n(ng),nnew)=max(cmu_fac3*0.5_r8* &
860 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
861 & (svstr(i,j)+svstr(i,j+1))**2), &
862 & gls_kmin(ng))
863# endif
864 tke(i,j,0,nnew)=max(cmu_fac3*0.5_r8* &
865 & sqrt((bustr(i,j)+bustr(i+1,j))**2+ &
866 & (bvstr(i,j)+bvstr(i,j+1))**2), &
867 & gls_kmin(ng))
868# if defined CHARNOK
869 zos_eff(i)=max(charnok_alpha(ng)/g*0.5_r8* &
870 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
871 & (svstr(i,j)+svstr(i,j+1))**2), &
872 & zos_min)
873# elif defined ZOS_HSIG
874 zos_eff(i)=max(zos_hsig_alpha(ng)*hwave(i,j), zos_min)
875# else
876 zos_eff(i)=zos_min
877# endif
878 gls(i,j,n(ng),nnew)=max(gls_cmu0(ng)**gls_p(ng)* &
879 & tke(i,j,n(ng),nnew)**gls_m(ng)* &
880 & (l_sft*zos_eff(i))**gls_n(ng), &
881 & gls_pmin(ng))
882 cff=gls_fac4*(vonkar*zob_min(i,j))**(gls_n(ng))
883 gls(i,j,0,nnew)=max(cff*tke(i,j,0,nnew)**(gls_m(ng)), &
884 & gls_pmin(ng))
885 END DO
886!
887! Solve tri-diagonal system for turbulent kinetic energy.
888!
889 DO i=istr,iend
890# if defined CRAIG_BANNER
891 tke_fluxt(i)=dt(ng)*crgban_cw(ng)* &
892 & (0.50_r8* &
893 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
894 & (svstr(i,j)+svstr(i,j+1))**2))**1.5_r8
895# elif defined TKE_WAVEDISS
896 tke_fluxt(i)=dt(ng)*sz_alpha(ng)*(dissip_break(i,j)+ &
897 & dissip_wcap(i,j))
898# else
899 tke_fluxt(i)=0.0_r8
900# endif
901 tke_fluxb(i)=0.0_r8
902!
903 cff=1.0_r8/bck(i,n(ng)-1)
904 cf(i,n(ng)-1)=cff*fck(i,n(ng)-1)
905 tke(i,j,n(ng)-1,nnew)=cff*(tke(i,j,n(ng)-1,nnew)+tke_fluxt(i))
906 END DO
907 DO i=istr,iend
908 DO k=n(ng)-2,1,-1
909 cff=1.0_r8/(bck(i,k)-cf(i,k+1)*fck(i,k+1))
910 cf(i,k)=cff*fck(i,k)
911 tke(i,j,k,nnew)=cff*(tke(i,j,k,nnew)- &
912 & fck(i,k+1)*tke(i,j,k+1,nnew))
913 END DO
914 tke(i,j,1,nnew)=tke(i,j,1,nnew)-cff*tke_fluxb(i)
915 END DO
916 DO k=2,n(ng)-1
917 DO i=istr,iend
918 tke(i,j,k,nnew)=tke(i,j,k,nnew)-cf(i,k)*tke(i,j,k-1,nnew)
919 END DO
920 END DO
921!
922! Solve tri-diagonal system for generic statistical field.
923!
924 DO i=istr,iend
925 cff=0.5_r8*(tke(i,j,n(ng),nnew)+tke(i,j,n(ng)-1,nnew))
926 gls_fluxt(i)=dt(ng)*gls_fac3*cff**gls_m(ng)* &
927 & l_sft**(gls_n(ng))* &
928 & (zos_eff(i)+0.5_r8*hz(i,j,n(ng)))** &
929 & (gls_n(ng)-1.0_r8)* &
930 & 0.5_r8*(akp(i,j,n(ng))+akp(i,j,n(ng)-1))
931# ifdef CRAIG_BANNER
932 gls_fluxt(i)=gls_fluxt(i)-dt(ng)* &
933 & gls_m(ng)*(gls_cmu0(ng)**gls_p(ng))* &
934 & cff**(gls_m(ng)-1.0_r8)* &
935 & ((zos_eff(i)+0.5_r8*hz(i,j,n(ng)))*l_sft)** &
936 & gls_n(ng)* &
937 & gls_sigk(ng)*ogls_sigp*crgban_cw(ng)* &
938 & (0.5_r8* &
939 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
940 & (svstr(i,j)+svstr(i,j+1))**2))**1.5_r8
941# elif defined TKE_WAVEDISS
942 gls_fluxt(i)=gls_fluxt(i)-dt(ng)* &
943 & gls_m(ng)*(gls_cmu0(ng)**gls_p(ng))* &
944 & cff**(gls_m(ng)-1.0_r8)* &
945 & ((zos_eff(i)+0.5_r8*hz(i,j,n(ng)))*l_sft)** &
946 & gls_n(ng)* &
947 & gls_sigk(ng)*ogls_sigp*sz_alpha(ng)* &
948 & (dissip_break(i,j)+dissip_wcap(i,j))
949# endif
950 cff=0.5_r8*(tke(i,j,0,nnew)+tke(i,j,1,nnew))
951 gls_fluxb(i)=dt(ng)*gls_fac2*(cff**gls_m(ng))* &
952 & (0.5_r8*hz(i,j,1)+zob_min(i,j))** &
953 & (gls_n(ng)-1.0_r8)* &
954 & 0.5_r8*(akp(i,j,0)+akp(i,j,1))
955!
956 cff=1.0_r8/bcp(i,n(ng)-1)
957 cf(i,n(ng)-1)=cff*fcp(i,n(ng)-1)
958 gls(i,j,n(ng)-1,nnew)=cff*(gls(i,j,n(ng)-1,nnew)-gls_fluxt(i))
959 END DO
960 DO i=istr,iend
961 DO k=n(ng)-2,1,-1
962 cff=1.0_r8/(bcp(i,k)-cf(i,k+1)*fcp(i,k+1))
963 cf(i,k)=cff*fcp(i,k)
964 gls(i,j,k,nnew)=cff*(gls(i,j,k,nnew)- &
965 & fcp(i,k+1)*gls(i,j,k+1,nnew))
966 END DO
967 gls(i,j,1,nnew)=gls(i,j,1,nnew)-cff*gls_fluxb(i)
968!! gls(i,j,1,nnew)=MAX(gls(i,j,1,nnew), gls_Pmin(ng))
969 END DO
970 DO k=2,n(ng)-1
971 DO i=istr,iend
972 gls(i,j,k,nnew)=gls(i,j,k,nnew)-cf(i,k)*gls(i,j,k-1,nnew)
973!! gls(i,j,k,nnew)=MAX(gls(i,j,k,nnew), gls_Pmin(ng))
974 END DO
975 END DO
976!
977!---------------------------------------------------------------------
978! Compute vertical mixing coefficients (m2/s).
979!---------------------------------------------------------------------
980!
981 DO i=istr,iend
982 DO k=1,n(ng)-1
983!
984! Compute turbulent length scale (m).
985!
986 tke(i,j,k,nnew)=max(tke(i,j,k,nnew),gls_kmin(ng))
987 gls(i,j,k,nnew)=max(gls(i,j,k,nnew),gls_pmin(ng))
988 IF (gls_n(ng).ge.0.0_r8) THEN
989 gls(i,j,k,nnew)=min(gls(i,j,k,nnew),gls_fac5* &
990 & tke(i,j,k,nnew)**(tke_exp4)* &
991 & (sqrt(max(0.0_r8, &
992 & buoy2(i,j,k)))+eps)** &
993 & (-gls_n(ng)))
994 ELSE
995 gls(i,j,k,nnew)=max(gls(i,j,k,nnew),gls_fac5* &
996 & tke(i,j,k,nnew)**(tke_exp4)* &
997 & (sqrt(max(0.0_r8, &
998 & buoy2(i,j,k)))+eps)** &
999 & (-gls_n(ng)))
1000 END IF
1001 ls_unlmt=max(eps, &
1002 & gls(i,j,k,nnew)**( gls_exp1)*cmu_fac1* &
1003 & tke(i,j,k,nnew)**(-tke_exp1))
1004 IF (buoy2(i,j,k).gt.0.0_r8) THEN
1005 ls_lmt=min(ls_unlmt, &
1006 & sqrt(0.56_r8*tke(i,j,k,nnew)/ &
1007 & (max(0.0_r8,buoy2(i,j,k))+eps)))
1008 ELSE
1009 ls_lmt=ls_unlmt
1010 END IF
1011!
1012! Recompute gls based on limited length scale
1013!
1014 gls(i,j,k,nnew)=max(gls_cmu0(ng)**gls_p(ng)* &
1015 & tke(i,j,k,nnew)**gls_m(ng)* &
1016 & ls_lmt**gls_n(ng), gls_pmin(ng))
1017!
1018! Compute nondimensional stability functions for tracers (Sh) and
1019! momentum (Sm).
1020!
1021 gh=min(gls_gh0,-buoy2(i,j,k)*ls_lmt*ls_lmt/ &
1022 & (2.0_r8*tke(i,j,k,nnew)))
1023 gh=min(gh,gh-(gh-gls_ghcri)**2/ &
1024 & (gh+gls_gh0-2.0_r8*gls_ghcri))
1025 gh=max(gh,gls_ghmin)
1026# if defined CANUTO_A || defined CANUTO_B
1027!
1028! Compute shear number.
1029!
1030 gm=(gls_b0/gls_fac6-gls_b1*gh+gls_b3*gls_fac6*(gh**2))/ &
1031 & (gls_b2-gls_b4*gls_fac6*gh)
1032 gm=min(gm,shear2(i,j,k)*ls_lmt*ls_lmt/ &
1033 & (2.0_r8*tke(i,j,k,nnew)))
1034!! Gm=MIN(Gm,(gls_s1*gls_fac6*Gh-gls_s0)/(gls_s2*gls_fac6))
1035!
1036! Compute stability functions
1037!
1038 cff=gls_b0-gls_b1*gls_fac6*gh+gls_b2*gls_fac6*gm+ &
1039 & gls_b3*gls_fac6**2*gh**2-gls_b4*gls_fac6**2*gh*gm+ &
1040 & gls_b5*gls_fac6**2*gm*gm
1041 sm=(gls_s0-gls_s1*gls_fac6*gh+gls_s2*gls_fac6*gm)/cff
1042 sh=(gls_s4-gls_s5*gls_fac6*gh+gls_s6*gls_fac6*gm)/cff
1043 sm=max(sm,0.0_r8)
1044 sh=max(sh,0.0_r8)
1045!
1046! Relate Canuto stability to ROMS notation
1047!
1048 sm=sm*sqrt2/gls_cmu0(ng)**3
1049 sh=sh*sqrt2/gls_cmu0(ng)**3
1050# elif defined KANTHA_CLAYSON
1051 cff=1.0_r8-my_sh2*gh
1052 sh=my_sh1/cff
1053 sm=(my_b1pm1o3+my_sm4*sh*gh)/(1.0_r8-my_sm2*gh)
1054# else /* Galperin */
1055 cff=1.0_r8-my_sh2*gh
1056 sh=my_sh1/cff
1057 sm=(my_sm3+sh*gh*my_sm4)/(1.0_r8-my_sm2*gh)
1058# endif
1059!
1060! Compute vertical mixing (m2/s) coefficients of momentum and
1061! tracers. Average ql over the two timesteps rather than using
1062! the new Lscale and just averaging tke.
1063!
1064 ql=sqrt2*0.5_r8*(ls_lmt*sqrt(tke(i,j,k,nnew))+ &
1065 & lscale(i,j,k)*sqrt(tke(i,j,k,nstp)))
1066 akv(i,j,k)=akv_bak(ng)+sm*ql
1067 DO itrc=1,nat
1068 akt(i,j,k,itrc)=akt_bak(itrc,ng)+sh*ql
1069 END DO
1070!
1071! Compute vertical mixing (m2/s) coefficents of turbulent kinetic
1072! energy and generic statistical field.
1073!
1074 akk(i,j,k)=akk_bak(ng)+ &
1075 & sm*ql/gls_sigk(ng)
1076
1077# if defined CRAIG_BANNER || defined TKE_WAVEDISS
1078!
1079! If wave breaking, modify surface boundary condition for
1080! gls diffusivity Schmidt number.
1081!
1082 pprod=gls_c1(ng)*shear2(i,j,k)*akv(i,j,k)
1083 cff=cmu_fac2*tke(i,j,k,nnew)**(1.5_r8+tke_exp1)* &
1084 & gls(i,j,k,nnew)**(-1.0_r8/gls_n(ng))
1085 cff2=min(pprod/cff, 1.0_r8)
1086 sig_eff=cff2*gls_sigp(ng)+(1.0_r8-cff2)*gls_sigp_cb
1087 akp(i,j,k)=akp_bak(ng)+sm*ql/sig_eff
1088# else
1089 akp(i,j,k)=akp_bak(ng)+sm*ql*ogls_sigp
1090# endif
1091!
1092! Save limited length scale.
1093!
1094 lscale(i,j,k)=ls_lmt
1095 END DO
1096!
1097! Compute vertical mixing coefficients at the surface and bottom.
1098!
1099!
1100 akv(i,j,n(ng))=akv_bak(ng)+l_sft*zos_eff(i)*gls_cmu0(ng)* &
1101 & sqrt(tke(i,j,n(ng),nnew))
1102 akv(i,j,0)=akv_bak(ng)+vonkar*zob_min(i,j)*gls_cmu0(ng)* &
1103 & sqrt(tke(i,j,0,nnew))
1104!
1105 akk(i,j,n(ng))=akk_bak(ng)+akv(i,j,n(ng))/gls_sigk(ng)
1106 akk(i,j,0)=akk_bak(ng)+akv(i,j,0)/gls_sigk(ng)
1107 akp(i,j,n(ng))=akp_bak(ng)+akv(i,j,n(ng))*ogls_sigp
1108 akp(i,j,0)=akp_bak(ng)+akv(i,j,0)/gls_sigp(ng)
1109!
1110 DO itrc=1,nat
1111 akt(i,j,n(ng),itrc)=akt_bak(itrc,ng)
1112 akt(i,j,0,itrc)=akt_bak(itrc,ng)
1113 END DO
1114
1115# if defined LIMIT_VDIFF || defined LIMIT_VVISC
1116!
1117! Limit vertical mixing coefficients with the upper threshold value.
1118! This is an engineering fix but it can be based on the fact that
1119! vertical mixing in the ocean from indirect observations are not
1120! higher than the threshold value.
1121!
1122 DO k=0,n(ng)
1123# ifdef LIMIT_VDIFF
1124 DO itrc=1,nat
1125 akt(i,j,k,itrc)=min(akt_limit(itrc,ng), akt(i,j,k,itrc))
1126 END DO
1127# endif
1128# ifdef LIMIT_VVISC
1129 akv(i,j,k)=min(akv_limit(ng), akv(i,j,k))
1130# endif
1131 END DO
1132# endif
1133 END DO
1134 END DO
1135!
1136!-----------------------------------------------------------------------
1137! Set lateral boundary conditions.
1138!-----------------------------------------------------------------------
1139!
1140 DO k=0,n(ng)
1141 IF (domain(ng)%Western_Edge(tile)) THEN
1142 DO j=jstr,jend
1143 DO itrc=1,nat
1144 akt(istr-1,j,k,itrc)=akt(istr,j,k,itrc)
1145 END DO
1146 akv(istr-1,j,k)=akv(istr,j,k)
1147 END DO
1148 END IF
1149 IF (domain(ng)%Eastern_Edge(tile)) THEN
1150 DO j=jstr,jend
1151 DO itrc=1,nat
1152 akt(iend+1,j,k,itrc)=akt(iend,j,k,itrc)
1153 END DO
1154 akv(iend+1,j,k)=akv(iend,j,k)
1155 END DO
1156
1157 END IF
1158 IF (domain(ng)%Southern_Edge(tile)) THEN
1159 DO i=istr,iend
1160 DO itrc=1,nat
1161 akt(i,jstr-1,k,itrc)=akt(i,jstr,k,itrc)
1162 END DO
1163 akv(i,jstr-1,k)=akv(i,jstr,k)
1164 END DO
1165 END IF
1166 IF (domain(ng)%Northern_Edge(tile)) THEN
1167 DO i=istr,iend
1168 DO itrc=1,nat
1169 akt(i,jend+1,k,itrc)=akt(i,jend,k,itrc)
1170 END DO
1171 akv(i,jend+1,k)=akv(i,jend,k)
1172 END DO
1173 END IF
1174 IF (domain(ng)%SouthWest_Corner(tile)) THEN
1175 DO itrc=1,nat
1176 akt(istr-1,jstr-1,k,itrc)=0.5_r8* &
1177 & (akt(istr ,jstr-1,k,itrc)+ &
1178 & akt(istr-1,jstr ,k,itrc))
1179 END DO
1180 akv(istr-1,jstr-1,k)=0.5_r8* &
1181 & (akv(istr ,jstr-1,k)+ &
1182 & akv(istr-1,jstr ,k))
1183 END IF
1184 IF (domain(ng)%SouthEast_Corner(tile)) THEN
1185 DO itrc=1,nat
1186 akt(iend+1,jstr-1,k,itrc)=0.5_r8* &
1187 & (akt(iend ,jstr-1,k,itrc)+ &
1188 & akt(iend+1,jstr ,k,itrc))
1189 END DO
1190 akv(iend+1,jstr-1,k)=0.5_r8* &
1191 & (akv(iend ,jstr-1,k)+ &
1192 & akv(iend+1,jstr ,k))
1193 END IF
1194 IF (domain(ng)%NorthWest_Corner(tile)) THEN
1195 DO itrc=1,nat
1196 akt(istr-1,jend+1,k,itrc)=0.5_r8* &
1197 & (akt(istr ,jend+1,k,itrc)+ &
1198 & akt(istr-1,jend ,k,itrc))
1199 END DO
1200 akv(istr-1,jend+1,k)=0.5_r8* &
1201 & (akv(istr ,jend+1,k)+ &
1202 & akv(istr-1,jend ,k))
1203 END IF
1204 IF (domain(ng)%NorthEast_Corner(tile)) THEN
1205 DO itrc=1,nat
1206 akt(iend+1,jend+1,k,itrc)=0.5_r8* &
1207 & (akt(iend ,jend+1,k,itrc)+ &
1208 & akt(iend+1,jend ,k,itrc))
1209 END DO
1210 akv(iend+1,jend+1,k)=0.5_r8* &
1211 & (akv(iend ,jend+1,k)+ &
1212 & akv(iend+1,jend ,k))
1213 END IF
1214 END DO
1215!
1216 CALL tkebc_tile (ng, tile, &
1217 & lbi, ubi, lbj, ubj, n(ng), &
1218 & imins, imaxs, jmins, jmaxs, &
1219 & nnew, nstp, &
1220 & gls, tke)
1221
1222 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
1223 CALL exchange_w3d_tile (ng, tile, &
1224 & lbi, ubi, lbj, ubj, 0, n(ng), &
1225 & tke(:,:,:,nnew))
1226 CALL exchange_w3d_tile (ng, tile, &
1227 & lbi, ubi, lbj, ubj, 0, n(ng), &
1228 & gls(:,:,:,nnew))
1229 CALL exchange_w3d_tile (ng, tile, &
1230 & lbi, ubi, lbj, ubj, 0, n(ng), &
1231 & akv)
1232 DO itrc=1,nat
1233 CALL exchange_w3d_tile (ng, tile, &
1234 & lbi, ubi, lbj, ubj, 0, n(ng), &
1235 & akt(:,:,:,itrc))
1236 END DO
1237 END IF
1238
1239# ifdef DISTRIBUTE
1240 CALL mp_exchange3d (ng, tile, inlm, 3, &
1241 & lbi, ubi, lbj, ubj, 0, n(ng), &
1242 & nghostpoints, &
1243 & ewperiodic(ng), nsperiodic(ng), &
1244 & tke(:,:,:,nnew), &
1245 & gls(:,:,:,nnew), &
1246 & akv)
1247 CALL mp_exchange4d (ng, tile, inlm, 1, &
1248 & lbi, ubi, lbj, ubj, 0, n(ng), 1, nat, &
1249 & nghostpoints, &
1250 & ewperiodic(ng), nsperiodic(ng), &
1251 & akt)
1252# endif
1253!
1254 RETURN
1255 END SUBROUTINE gls_corstep_tile
1256#endif
1257 END MODULE gls_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
gls_corstep_mod
Definition gls_corstep.F:2
gls_corstep_mod::gls_corstep_tile
subroutine gls_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, zobot, u, v, w, bustr, bvstr, sustr, svstr, hwave, dissip_break, dissip_wcap, akt, akv, bvf, akk, akp, lscale, gls, tke)
Definition gls_corstep.F:132
gls_corstep_mod::gls_corstep
subroutine, public gls_corstep(ng, tile)
Definition gls_corstep.F:38
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_b1pm1o3
real(r8) my_b1pm1o3
Definition mod_scalars.F:1862
mod_scalars::gls_b5
real(r8) gls_b5
Definition mod_scalars.F:1766
mod_scalars::gls_p
real(r8), dimension(:), allocatable gls_p
Definition mod_scalars.F:1741
mod_scalars::gls_n
real(r8), dimension(:), allocatable gls_n
Definition mod_scalars.F:1740
mod_scalars::gls_b3
real(r8) gls_b3
Definition mod_scalars.F:1764
mod_scalars::gls_m
real(r8), dimension(:), allocatable gls_m
Definition mod_scalars.F:1739
mod_scalars::gls_s2
real(r8) gls_s2
Definition mod_scalars.F:1756
mod_scalars::zos_hsig_alpha
real(r8), dimension(:), allocatable zos_hsig_alpha
Definition mod_scalars.F:1805
mod_scalars::gls_s4
real(r8) gls_s4
Definition mod_scalars.F:1758
mod_scalars::vonkar
real(dp) vonkar
Definition mod_scalars.F:469
mod_scalars::gls_ghcri
real(r8), parameter gls_ghcri
Definition mod_scalars.F:1770
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::gls_b0
real(r8) gls_b0
Definition mod_scalars.F:1761
mod_scalars::sz_alpha
real(r8), dimension(:), allocatable sz_alpha
Definition mod_scalars.F:1806
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::gls_b2
real(r8) gls_b2
Definition mod_scalars.F:1763
mod_scalars::gls_c2
real(r8), dimension(:), allocatable gls_c2
Definition mod_scalars.F:1747
mod_scalars::charnok_alpha
real(r8), dimension(:), allocatable charnok_alpha
Definition mod_scalars.F:1804
mod_scalars::gls_b4
real(r8) gls_b4
Definition mod_scalars.F:1765
mod_scalars::zos
real(r8), dimension(:), allocatable zos
Definition mod_scalars.F:777
mod_scalars::akk_bak
real(r8), dimension(:), allocatable akk_bak
Definition mod_scalars.F:889
mod_scalars::gls_cmu0
real(r8), dimension(:), allocatable gls_cmu0
Definition mod_scalars.F:1744
mod_scalars::gls_ghmin
real(r8), parameter gls_ghmin
Definition mod_scalars.F:1794
mod_scalars::gls_sigk
real(r8), dimension(:), allocatable gls_sigk
Definition mod_scalars.F:1742
mod_scalars::my_sm3
real(r8) my_sm3
Definition mod_scalars.F:1868
mod_scalars::crgban_cw
real(r8), dimension(:), allocatable crgban_cw
Definition mod_scalars.F:1807
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::gls_pmin
real(r8), dimension(:), allocatable gls_pmin
Definition mod_scalars.F:1751
mod_scalars::gls_e2
real(r8), parameter gls_e2
Definition mod_scalars.F:1795
mod_scalars::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::gls_s1
real(r8) gls_s1
Definition mod_scalars.F:1755
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::akp_bak
real(r8), dimension(:), allocatable akp_bak
Definition mod_scalars.F:890
mod_scalars::akt_limit
real(r8), dimension(:,:), allocatable akt_limit
Definition mod_scalars.F:908
mod_scalars::gls_c1
real(r8), dimension(:), allocatable gls_c1
Definition mod_scalars.F:1746
mod_scalars::akv_limit
real(r8), dimension(:), allocatable akv_limit
Definition mod_scalars.F:892
mod_scalars::gls_s6
real(r8) gls_s6
Definition mod_scalars.F:1760
mod_scalars::gls_gh0
real(r8), parameter gls_gh0
Definition mod_scalars.F:1769
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::gls_sigp
real(r8), dimension(:), allocatable gls_sigp
Definition mod_scalars.F:1743
mod_scalars::gls_s5
real(r8) gls_s5
Definition mod_scalars.F:1759
mod_scalars::gls_s0
real(r8) gls_s0
Definition mod_scalars.F:1754
mod_scalars::gls_kmin
real(r8), dimension(:), allocatable gls_kmin
Definition mod_scalars.F:1750
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::gls_b1
real(r8) gls_b1
Definition mod_scalars.F:1762
mod_scalars::gls_c3m
real(r8), dimension(:), allocatable gls_c3m
Definition mod_scalars.F:1748
mod_scalars::gls_c3p
real(r8), dimension(:), allocatable gls_c3p
Definition mod_scalars.F:1749
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
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