ROMS
Loading...
Searching...
No Matches
lmd_bkpp.F
Go to the documentation of this file.
1#include "cppdefs.h"
2#define SASHA
3 MODULE lmd_bkpp_mod
4#if defined NONLINEAR && defined LMD_BKPP && defined SOLVE3D
5!
6!git $Id$
7!================================================== Hernan G. Arango ===
8! Copyright (c) 2002-2025 The ROMS Group Scott M. Durski !
9! Licensed under a MIT/X style license !
10! See License_ROMS.md !
11!=======================================================================
12! !
13! This routine determines the depth of bottom oceanic boundary !
14! layer, hbbl, as the deepest depth where the bulk Richardson !
15! number is equal to the critical value, Ric. !
16! !
17! Then, it computes the vertical mixing coefficients within the !
18! boundary layer. They depend on surface forcing and the magnitude !
19! and gradient of interior mixing below the boundary layer. The !
20! ocean interior is allowed to force the boundary layer through a !
21! dependence of the nondimensional vertical shape function G(sigma) !
22! and its vertical derivative at sigma=1 on the interior mixing !
23! coefficients, and it vertical derivative at d=hsbl. The boundary !
24! layer mixing coefficients are computed by matching these values. !
25! !
26! Reference: !
27! !
28! Large, W.G., J.C. McWilliams, and S.C. Doney, 1994: A Review !
29! and model with a nonlocal boundary layer parameterization, !
30! Reviews of Geophysics, 32,363-403. !
31! !
32! This routine was adapted from Bill Large 1995 code. !
33! !
34!=======================================================================
35!
36 implicit none
37
38 PRIVATE
39 PUBLIC :: lmd_bkpp
40
41 CONTAINS
42!
43!***********************************************************************
44 SUBROUTINE lmd_bkpp (ng, tile)
45!***********************************************************************
46!
47 USE mod_param
48 USE mod_forces
49 USE mod_grid
50 USE mod_mixing
51 USE mod_ocean
52 USE mod_stepping
53!
54! Imported variable declarations.
55!
56 integer, intent(in) :: ng, tile
57!
58! Local variable declarations.
59!
60# include "tile.h"
61!
62 CALL lmd_bkpp_tile (ng, tile, &
63 & lbi, ubi, lbj, ubj, &
64 & imins, imaxs, jmins, jmaxs, &
65 & nstp(ng), &
66# ifdef MASKING
67 & grid(ng) % rmask, &
68# endif
69 & grid(ng) % f, &
70 & grid(ng) % h, &
71 & grid(ng) % Hz, &
72 & grid(ng) % z_r, &
73 & grid(ng) % z_w, &
74 & ocean(ng) % u, &
75 & ocean(ng) % v, &
76 & ocean(ng) % pden, &
77 & forces(ng) % srflx, &
78 & forces(ng) % btflx, &
79 & forces(ng) % bustr, &
80 & forces(ng) % bvstr, &
81 & mixing(ng) % alpha, &
82# ifdef SALINITY
83 & mixing(ng) % beta, &
84# endif
85 & mixing(ng) % bvf, &
86 & mixing(ng) % ksbl, &
87 & mixing(ng) % Akt, &
88 & mixing(ng) % Akv, &
89 & mixing(ng) % kbbl, &
90 & mixing(ng) % hbbl)
91 RETURN
92 END SUBROUTINE lmd_bkpp
93!
94!***********************************************************************
95 SUBROUTINE lmd_bkpp_tile (ng, tile, &
96 & LBi, UBi, LBj, UBj, &
97 & IminS, ImaxS, JminS, JmaxS, &
98 & nstp, &
99# ifdef MASKING
100 & rmask, &
101# endif
102 & f, h, Hz, z_r, z_w, &
103 & u, v, pden, &
104 & srflx, btflx, bustr, bvstr, &
105 & alpha, &
106# ifdef SALINITY
107 & beta, &
108# endif
109 & bvf, &
110 & ksbl, Akt, Akv, kbbl, hbbl)
111!***********************************************************************
112!
113 USE mod_param
114 USE mod_scalars
115!
116 USE bc_2d_mod, ONLY : bc_r2d_tile
117# ifdef DISTRIBUTE
118 USE mp_exchange_mod, ONLY : mp_exchange2d
119# endif
120# ifdef LMD_SHAPIRO
121 USE shapiro_mod
122# endif
123!
124! Imported variable declarations.
125!
126 integer, intent(in) :: ng, tile
127 integer, intent(in) :: LBi, UBi, LBj, UBj
128 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
129 integer, Intent(in) :: nstp
130!
131# ifdef ASSUMED_SHAPE
132# ifdef MASKING
133 real(r8), intent(in) :: rmask(LBi:,LBj:)
134# endif
135 real(r8), intent(in) :: f(LBi:,LBj:)
136 real(r8), intent(in) :: h(LBi:,LBj:)
137 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
138 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
139 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
140 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
141 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
142 real(r8), intent(in) :: pden(LBi:,LBj:,:)
143 real(r8), intent(in) :: srflx(LBi:,LBj:)
144 real(r8), intent(in) :: btflx(LBi:,LBj:,:)
145 real(r8), intent(in) :: bustr(LBi:,LBj:)
146 real(r8), intent(in) :: bvstr(LBi:,LBj:)
147 real(r8), intent(in) :: alpha(LBi:,LBj:)
148# ifdef SALINITY
149 real(r8), intent(in) :: beta(LBi:,LBj:)
150# endif
151 real(r8), intent(in) :: bvf(LBi:,LBj:,0:)
152 integer, intent(in) :: ksbl(LBi:,LBj:)
153 real(r8), intent(inout) :: Akt(LBi:,LBj:,0:,:)
154 real(r8), intent(inout) :: Akv(LBi:,LBj:,0:)
155 real(r8), intent(inout) :: hbbl(LBi:,LBj:)
156 integer, intent(out) :: kbbl(LBi:,LBj:)
157# else
158# ifdef MASKING
159 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
160# endif
161 real(r8), intent(in) :: f(LBi:UBi,LBj:UBj)
162 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
163 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
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),3)
167 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),3)
168 real(r8), intent(in) :: pden(LBi:UBi,LBj:UBj,N(ng))
169 real(r8), intent(in) :: srflx(LBi:UBi,LBj:UBj)
170 real(r8), intent(in) :: btflx(LBi:UBi,LBj:UBj,NT(ng))
171 real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj)
172 real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
173 real(r8), intent(in) :: alpha(LBi:UBi,LBj:UBj)
174# ifdef SALINITY
175 real(r8), intent(in) :: beta(LBi:UBi,LBj:UBj)
176# endif
177 real(r8), intent(in) :: bvf(LBi:UBi,LBj:UBj,0:N(ng))
178 integer, intent(in) :: ksbl(LBi:UBi,LBj:UBj)
179 real(r8), intent(inout) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
180 real(r8), intent(inout) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
181 real(r8), intent(inout) :: hbbl(LBi:UBi,LBj:UBj)
182 integer, intent(out) :: kbbl(LBi:UBi,LBj:UBj)
183# endif
184!
185! Local variable declarations.
186!
187 integer :: i, itrc, j, k
188
189 real(r8), parameter :: eps = 1.0e-10_r8
190 real(r8), parameter :: r3 = 1.0_r8/3.0_r8
191 real(r8), parameter :: small = 1.0e-20_r8
192
193 real(r8) :: Gm, Gt, Gs, K_bl, Ribot, Ritop, Rk
194 real(r8) :: Uk, Ustar3, Vk, Vtc
195 real(r8) :: a1, a2, a3, cff, cff1, cff2, cff_up, cff_dn
196 real(r8) :: depth, dK_bl, hekman, sigma, zbl
197 real(r8) :: zetahat, zetapar
198
199 real(r8), dimension (IminS:ImaxS) :: Rref
200 real(r8), dimension (IminS:ImaxS) :: Uref
201 real(r8), dimension (IminS:ImaxS) :: Vref
202
203 real(r8), dimension (IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: Bflux
204
205 real(r8), dimension (IminS:ImaxS,0:N(ng)) :: FC
206 real(r8), dimension (IminS:ImaxS,0:N(ng)) :: dR
207 real(r8), dimension (IminS:ImaxS,0:N(ng)) :: dU
208 real(r8), dimension (IminS:ImaxS,0:N(ng)) :: dV
209
210 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: Bo
211 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: Bosol
212 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: Bfbot
213 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: Gm1
214 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: Gt1
215 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: Gs1
216 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: Ustar
217 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: bl_dpth
218 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: dGm1dS
219 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: dGt1dS
220 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: dGs1dS
221 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: f1
222 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: swdk
223 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: wm
224 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: ws
225 real(r8), dimension (IminS:ImaxS,JminS:JmaxS) :: zgrid
226
227# include "set_bounds.h"
228!
229!-----------------------------------------------------------------------
230! Initialize relevant parameters.
231!-----------------------------------------------------------------------
232!
233 vtc=lmd_cv*sqrt(-lmd_betat)/(sqrt(lmd_cs*lmd_epsilon)* &
234 & lmd_ric*vonkar*vonkar)
235!
236!-----------------------------------------------------------------------
237! Get approximation of bottom layer depth using "lmd_eps" and
238! boundary layer depth from previous time step.
239!-----------------------------------------------------------------------
240!
241 DO j=jstr,jend
242 DO i=istr,iend
243 bl_dpth(i,j)=lmd_epsilon*(hbbl(i,j)-z_w(i,j,0))
244 END DO
245 END DO
246!
247!-----------------------------------------------------------------------
248! Compute turbulent friction velocity (m/s) "Ustar" from bottom
249! stress at RHO-points.
250!-----------------------------------------------------------------------
251!
252 DO j=jstr,jend
253 DO i=istr,iend
254 ustar(i,j)=sqrt(sqrt((0.5_r8*(bustr(i,j)+bustr(i+1,j)))**2+ &
255 & (0.5_r8*(bvstr(i,j)+bvstr(i,j+1)))**2))
256# ifdef MASKING
257 ustar(i,j)=ustar(i,j)*rmask(i,j)
258# endif
259 END DO
260 END DO
261!
262!-----------------------------------------------------------------------
263! Compute bottom turbulent buoyancy forcing "Bo" (m2/s3). Compute
264! surface radiative buoyancy forcing "Bosol" (m2/s3) that can be
265! used in shallow areas when it can penetrate all the way to the
266! bottom.
267!-----------------------------------------------------------------------
268!
269 DO j=jstr,jend
270 DO i=istr,iend
271# ifdef SALINITY
272 bo(i,j)=g*(alpha(i,j)*btflx(i,j,itemp)- &
273 & beta(i,j)*btflx(i,j,isalt))
274# else
275 bo(i,j)=g*alpha(i,j)*btflx(i,j,itemp)
276# endif
277 bosol(i,j)=g*alpha(i,j)*srflx(i,j)
278 END DO
279 END DO
280!
281!-----------------------------------------------------------------------
282! Compute total buoyancy flux (m2/s3) at W-points.
283!-----------------------------------------------------------------------
284!
285 DO k=0,n(ng)
286 DO j=jstr,jend
287 DO i=istr,iend
288 zgrid(i,j)=z_w(i,j,n(ng))-z_w(i,j,k)
289 END DO
290 END DO
291 CALL lmd_swfrac_tile (ng, tile, &
292 & lbi, ubi, lbj, ubj, &
293 & imins, imaxs, jmins, jmaxs, &
294 & -1.0_r8, zgrid, swdk)
295 DO j=jstr,jend
296 DO i=istr,iend
297 bflux(i,j,k)=(bo(i,j)+bosol(i,j)*(1.0_r8-swdk(i,j)))
298# ifdef MASKING
299 bflux(i,j,k)=bflux(i,j,k)*rmask(i,j)
300# endif
301 END DO
302 END DO
303 END DO
304!
305!=======================================================================
306! Compute bulk Richardson number "Rib" and then find depth of the
307! oceanic bottom boundary layer "hbbl", such that Rib(hbbl)=Ric.
308!=======================================================================
309!
310 DO j=jstr,jend
311# ifdef RI_SPLINES
312!
313! Construct parabolic splines for vertical derivatives of potential
314! density and velocity components at W-points. FC is a scratch array.
315!
316 DO i=istr,iend
317 fc(i,0)=0.0_r8
318 dr(i,0)=0.0_r8
319 du(i,0)=0.0_r8
320 dv(i,0)=0.0_r8
321 END DO
322 DO k=1,n(ng)-1
323 DO i=istr,iend
324 cff=1.0_r8/(2.0_r8*hz(i,j,k+1)+ &
325 & hz(i,j,k)*(2.0_r8-fc(i,k-1)))
326 fc(i,k)=cff*hz(i,j,k+1)
327 dr(i,k)=cff*(6.0_r8*(pden(i,j,k+1)-pden(i,j,k))- &
328 & hz(i,j,k)*dr(i,k-1))
329 du(i,k)=cff*(3.0_r8*(u(i ,j,k+1,nstp)-u(i, j,k,nstp)+ &
330 & u(i+1,j,k+1,nstp)-u(i+1,j,k,nstp))- &
331 & hz(i,j,k)*du(i,k-1))
332 dv(i,k)=cff*(3.0_r8*(v(i,j ,k+1,nstp)-v(i,j ,k,nstp)+ &
333 & v(i,j+1,k+1,nstp)-v(i,j+1,k,nstp))- &
334 & hz(i,j,k)*dv(i,k-1))
335 END DO
336 END DO
337 DO i=istr,iend
338 dr(i,n(ng))=0.0_r8
339 du(i,n(ng))=0.0_r8
340 dv(i,n(ng))=0.0_r8
341 END DO
342 DO k=n(ng)-1,1,-1
343 DO i=istr,iend
344 dr(i,k)=dr(i,k)-fc(i,k)*dr(i,k+1)
345 du(i,k)=du(i,k)-fc(i,k)*du(i,k+1)
346 dv(i,k)=dv(i,k)-fc(i,k)*dv(i,k+1)
347 END DO
348 END DO
349# else
350!
351! Compute vertical derivatives of potential density and velocity
352! components at W-points.
353!
354 DO k=1,n(ng)-1
355 DO i=istr,iend
356 cff=1.0_r8/(z_r(i,j,k+1)-z_r(i,j,k))
357 dr(i,k)=cff*(pden(i,j,k+1)-pden(i,j,k))
358 cff=0.5_r8*cff
359 du(i,k)=cff*(u(i ,j,k+1,nstp)-u(i, j,k,nstp)+ &
360 & u(i+1,j,k+1,nstp)-u(i+1,j,k,nstp))
361 dv(i,k)=cff*(v(i,j ,k+1,nstp)-v(i,j ,k,nstp)+ &
362 & v(i,j+1,k+1,nstp)-v(i,j+1,k,nstp))
363 END DO
364 END DO
365 DO i=istr,iend
366 dr(i,0)=0.0_r8
367 dr(i,n(ng))=0.0_r8
368 du(i,0)=0.0_r8
369 du(i,n(ng))=0.0_r8
370 dv(i,0)=0.0_r8
371 dv(i,n(ng))=0.0_r8
372 END DO
373# endif
374!
375!-----------------------------------------------------------------------
376! Compute bulk Richardson number "Rib" and then find depth of oceanic
377! bottom boundary layer "hbbl".
378!
379! [Br - B(d)] * d
380! Rib(d) = ----------------------- ; Rib(hbbl)=Ric (1)
381! |Vr - V(d)|^2 + Vt(d)^2
382!
383! where "Br" and "Vr" are the bottom reference buoyancy and velocity
384! while "B(d)" and "V(d)" are the buoyancy and velocity at depth "d".
385!
386! In the code below, the criterion "Rib(hbbl)=Ric" is reformulated
387! as follows:
388!
389! Rib(d) Ritop(d)
390! ------ = --------------- = 1 (2)
391! Ric Ric * Ribot(d)
392!
393! where "Ritop" and "Ribot" are numerator and denominator in Eq. (1).
394! In its turn, Eq. (2) is rewritten in the following form:
395!
396! FC(d) = Ritop(d) - Ric * Ribot(d) = 0 (3)
397!
398! That is, the planetary boundary layer extends to the depth where
399! the critical function "FC(d)" changes its sign.
400!-----------------------------------------------------------------------
401!
402! Compute potential density and velocity components bottom reference
403! values.
404!
405 cff1=1.0_r8/3.0_r8
406 cff2=1.0_r8/6.0_r8
407 DO i=istr,iend
408 rref(i)=pden(i,j,1)- &
409 & hz(i,j,1)*(cff1*dr(i,0)+cff2*dr(i,1))
410 uref(i)=0.5_r8*(u(i,j,1,nstp)+u(i+1,j,1,nstp))- &
411 & hz(i,j,1)*(cff1*du(i,0)+cff2*du(i,1))
412 vref(i)=0.5_r8*(v(i,j,1,nstp)+v(i,j+1,1,nstp))- &
413 & hz(i,j,1)*(cff1*dv(i,0)+cff2*dv(i,1))
414 END DO
415!
416! Compute turbulent velocity scales for momentum (wm) and tracers (ws).
417! Then, compute critical function (FC) for bulk Richardson number.
418!
419 DO i=istr,iend
420 fc(i,0)=0.0_r8
421 DO k=1,n(ng)
422 depth=z_w(i,j,k)-z_w(i,j,0)
423 IF (bflux(i,j,k).lt.0.0_r8) THEN
424 sigma=min(bl_dpth(i,j),depth)
425 ELSE
426 sigma=depth
427 END IF
428 ustar3=ustar(i,j)*ustar(i,j)*ustar(i,j)
429 zetahat=vonkar*sigma*bflux(i,j,k)
430 zetapar=zetahat/(ustar3+small)
431 IF (zetahat.ge.0.0_r8) THEN ! stable
432 wm(i,j)=vonkar*ustar(i,j)/(1.0_r8+5.0_r8*zetapar)
433 ws(i,j)=wm(i,j)
434 ELSE ! unstable
435 IF (zetapar.gt.lmd_zetam) THEN
436 wm(i,j)=vonkar*ustar(i,j)* &
437 & (1.0_r8-16.0_r8*zetapar)**0.25_r8
438 ELSE
439 wm(i,j)=vonkar*(lmd_am*ustar3-lmd_cm*zetahat)**r3
440 END IF
441 IF (zetapar.gt.lmd_zetas) THEN
442 ws(i,j)=vonkar*ustar(i,j)* &
443 & (1.0_r8-16.0_r8*zetapar)**0.5_r8
444 ELSE
445 ws(i,j)=vonkar*(lmd_as*ustar3-lmd_cs*zetahat)**r3
446 END IF
447 END IF
448!
449 rk=pden(i,j,k)+ &
450 & hz(i,j,k)*(cff1*dr(i,k)+cff2*dr(i,k-1))
451 uk=0.5_r8*(u(i,j,k,nstp)+u(i+1,j,k,nstp))+ &
452 & hz(i,j,k)*(cff1*du(i,k)+cff2*du(i,k-1))
453 vk=0.5_r8*(v(i,j,k,nstp)+v(i,j+1,k,nstp))+ &
454 & hz(i,j,k)*(cff1*dv(i,k)+cff2*dv(i,k-1))
455!
456 ritop=-gorho0*(rk-rref(i))*depth
457 ribot=(uk-uref(i))**2+(vk-vref(i))**2+ &
458 & vtc*depth*ws(i,j)*sqrt(abs(bvf(i,j,k)))
459# ifdef SASHA
460 fc(i,k)=ritop-lmd_ric*ribot
461# else
462 fc(i,k)=ritop/(ribot+eps)
463# endif
464 END DO
465 END DO
466!
467! Linearly interpolate to find "hbbl" where Rib/Ric=1.
468!
469 DO i=istr,iend
470 kbbl(i,j)=n(ng)
471 hbbl(i,j)=z_w(i,j,n(ng))
472 END DO
473# ifdef SASHA
474 DO k=1,n(ng)-1
475 DO i=istr,iend
476 IF ((kbbl(i,j).eq.n(ng)).and.(fc(i,k).gt.0.0_r8)) THEN
477 hbbl(i,j)=(z_w(i,j,k)*fc(i,k-1)-z_w(i,j,k-1)*fc(i,k))/ &
478 & (fc(i,k-1)-fc(i,k))
479 kbbl(i,j)=k
480 END IF
481 END DO
482 END DO
483# else
484 DO k=1,n(ng)
485 DO i=istr,iend
486 IF ((kbbl(i,j).eq.n(ng)).and.((fc(i,k-1).lt.lmd_ric).and. &
487 & (fc(i,k ).ge.lmd_ric))) THEN
488 hbbl(i,j)=((lmd_ric-fc(i,k-1))*z_w(i,j,k )+ &
489 & (fc(i,k )-lmd_ric)*z_w(i,j,k-1))/ &
490 & (fc(i,k)-fc(i,k-1))
491 kbbl(i,j)=k
492 END IF
493 END DO
494 END DO
495# endif
496 END DO
497!
498! Compute total buoyancy flux at bottom boundary layer depth,
499! "Bfbot".
500!
501 DO j=jstr,jend
502 DO i=istr,iend
503 zgrid(i,j)=z_w(i,j,n(ng))-hbbl(i,j)
504# ifdef MASKING
505 zgrid(i,j)=zgrid(i,j)*rmask(i,j)
506# endif
507 END DO
508 END DO
509 CALL lmd_swfrac_tile (ng, tile, &
510 & lbi, ubi, lbj, ubj, &
511 & imins, imaxs, jmins, jmaxs, &
512 & -1.0_r8, zgrid, swdk)
513 DO j=jstr,jend
514 DO i=istr,iend
515 bfbot(i,j)=(bo(i,j)+bosol(i,j)*(1.0_r8-swdk(i,j)))
516# ifdef MASKING
517 bfbot(i,j)=bfbot(i,j)*rmask(i,j)
518# endif
519 END DO
520 END DO
521!
522! Under neutral and stable conditions, the depth of the bottom
523! boundary layer is required to be less than Ekman and Monin-Obukov
524! depths.
525!
526 DO j=jstr,jend
527 DO i=istr,iend
528 IF (ustar(i,j).ge.0.0_r8) THEN
529 hekman=lmd_cekman*ustar(i,j)/max(abs(f(i,j)),eps)-h(i,j)
530 hbbl(i,j)= min(hekman,hbbl(i,j))
531 END IF
532 hbbl(i,j)=min(hbbl(i,j),z_w(i,j,n(ng)))
533 hbbl(i,j)=max(hbbl(i,j),z_w(i,j,0))
534# ifdef MASKING
535 hbbl(i,j)=hbbl(i,j)*rmask(i,j)
536# endif
537 END DO
538 END DO
539# ifdef LMD_SHAPIRO
540!
541! Apply gradient or periodic boundary conditions.
542!
543 CALL bc_r2d_tile (ng, tile, &
544 & lbi, ubi, lbj, ubj, &
545 & hbbl)
546# ifdef DISTRIBUTE
547 CALL mp_exchange2d (ng, tile, inlm, 1, &
548 & lbi, ubi, lbj, ubj, &
549 & nghostpoints, &
550 & ewperiodic(ng), nsperiodic(ng), &
551 & hbbl)
552# endif
553!
554! Shapiro filter on boundary layer thickness
555!
556 CALL shapiro2d_tile (ng, tile, inlm, &
557 & lbi, ubi, lbj, ubj, &
558 & imins, imaxs, jmins, jmaxs, &
559# ifdef MASKING
560 & rmask, &
561# endif
562 & hbbl)
563!
564 DO j=jstr,jend
565 DO i=istr,iend
566 hbbl(i,j)=min(hbbl(i,j),z_w(i,j,n(ng)))
567 hbbl(i,j)=max(hbbl(i,j),z_w(i,j,0))
568# ifdef MASKING
569 hbbl(i,j)=hbbl(i,j)*rmask(i,j)
570# endif
571 END DO
572 END DO
573# endif
574!
575! Apply gradient or periodic boundary conditions.
576!
577 CALL bc_r2d_tile (ng, tile, &
578 & lbi, ubi, lbj, ubj, &
579 & hbbl)
580# ifdef DISTRIBUTE
581 CALL mp_exchange2d (ng, tile, inlm, 1, &
582 & lbi, ubi, lbj, ubj, &
583 & nghostpoints, &
584 & ewperiodic(ng), nsperiodic(ng), &
585 & hbbl)
586# endif
587!
588! Find new boundary layer index "kbbl".
589!
590 DO j=jstr,jend
591 DO i=istr,iend
592 kbbl(i,j)=n(ng)
593 DO k=1,n(ng)
594 IF ((kbbl(i,j).eq.n(ng)).and.(z_w(i,j,k).gt.hbbl(i,j))) THEN
595 kbbl(i,j)=k
596 END IF
597 END DO
598 END DO
599 END DO
600!
601!-----------------------------------------------------------------------
602! Compute total buoyancy flux at final bottom boundary layer depth.
603!-----------------------------------------------------------------------
604!
605 DO j=jstr,jend
606 DO i=istr,iend
607 zgrid(i,j)=z_w(i,j,n(ng))-hbbl(i,j)
608# ifdef MASKING
609 zgrid(i,j)=zgrid(i,j)*rmask(i,j)
610# endif
611 END DO
612 END DO
613 CALL lmd_swfrac_tile (ng, tile, &
614 & lbi, ubi, lbj, ubj, &
615 & imins, imaxs, jmins, jmaxs, &
616 & -1.0_r8, zgrid, swdk)
617 DO j=jstr,jend
618 DO i=istr,iend
619 bfbot(i,j)=(bo(i,j)+bosol(i,j)*(1.0_r8-swdk(i,j)))
620# ifdef MASKING
621 bfbot(i,j)=bfbot(i,j)*rmask(i,j)
622# endif
623 END DO
624 END DO
625!
626!=======================================================================
627! Compute vertical mixing coefficients within the planetary boundary
628! layer.
629!=======================================================================
630!
631! Compute tubulent velocity scales (wm,ws) at "hbbl".
632!
633 DO j=jstr,jend
634 DO i=istr,iend
635 bl_dpth(i,j)=lmd_epsilon*(hbbl(i,j)-z_w(i,j,0))
636 IF (bfbot(i,j).gt.0.0_r8) THEN
637 cff=1.0_r8
638 ELSE
639 cff=lmd_epsilon
640 END IF
641 sigma=cff*(hbbl(i,j)-z_w(i,j,0))
642 ustar3=ustar(i,j)*ustar(i,j)*ustar(i,j)
643 zetahat=vonkar*sigma*bfbot(i,j)
644 zetapar=zetahat/(ustar3+small)
645 IF (zetahat.ge.0.0_r8) THEN ! stable
646 wm(i,j)=vonkar*ustar(i,j)/(1.0_r8+5.0_r8*zetapar)
647 ws(i,j)=wm(i,j)
648 ELSE ! unstable
649 IF (zetapar.gt.lmd_zetam) THEN
650 wm(i,j)=vonkar*ustar(i,j)* &
651 & (1.0_r8-16.0_r8*zetapar)**0.25_r8
652 ELSE
653 wm(i,j)=vonkar*(lmd_am*ustar3-lmd_cm*zetahat)**r3
654 END IF
655 IF (zetapar.gt.lmd_zetas) THEN
656 ws(i,j)=vonkar*ustar(i,j)* &
657 & (1.0_r8-16.0_r8*zetapar)**0.5_r8
658 ELSE
659 ws(i,j)=vonkar*(lmd_as*ustar3-lmd_cs*zetahat)**r3
660 END IF
661 END IF
662 END DO
663 END DO
664!
665!-----------------------------------------------------------------------
666! Compute nondimensional shape function Gx(sigma) in terms of the
667! interior diffusivities at sigma=1 (Gm1, Gs1, Gt1) and its vertical
668! derivative evaluated "hbbl" via interpolation.
669!-----------------------------------------------------------------------
670!
671 DO j=jstr,jend
672 DO i=istr,iend
673 f1(i,j)=5.0_r8*max(0.0_r8,bfbot(i,j))*vonkar/ &
674 & (ustar(i,j)*ustar(i,j)*ustar(i,j)*ustar(i,j)+eps)
675 END DO
676 END DO
677!
678 DO j=jstr,jend
679 DO i=istr,iend
680 zbl=hbbl(i,j)-z_w(i,j,0)
681 k=kbbl(i,j)
682 cff=1.0_r8/(z_w(i,j,k)-z_w(i,j,k-1))
683 cff_dn=cff*(hbbl(i,j)-z_w(i,j,k-1))
684 cff_up=cff*(z_w(i,j,k)-hbbl(i,j))
685!
686! Compute nondimensional shape function for viscosity "Gm1" and its
687! vertical derivative "dGm1dS" evaluated at "hbbl".
688!
689 k_bl=cff_dn*akv(i,j,k)+cff_up*akv(i,j,k-1)
690 dk_bl=-cff*(akv(i,j,k)-akv(i,j,k-1))
691 gm1(i,j)=k_bl/(zbl*wm(i,j)+eps)
692# ifdef MASKING
693 gm1(i,j)=gm1(i,j)*rmask(i,j)
694# endif
695 dgm1ds(i,j)=min(0.0_r8, k_bl*f1(i,j)-dk_bl/(wm(i,j)+eps))
696!
697! Compute nondimensional shape function for diffusion of temperature
698! "Gt1" and its vertical derivative "dGt1dS" evaluated at "hbbl".
699!
700 k_bl=cff_dn*akt(i,j,k,itemp)+cff_up*akt(i,j,k-1,itemp)
701 dk_bl=-cff*(akt(i,j,k,itemp)-akt(i,j,k-1,itemp))
702 gt1(i,j)=k_bl/(zbl*ws(i,j)+eps)
703# ifdef MASKING
704 gt1(i,j)=gt1(i,j)*rmask(i,j)
705# endif
706 dgt1ds(i,j)=min(0.0_r8, k_bl*f1(i,j)-dk_bl/(ws(i,j)+eps))
707# ifdef SALINITY
708!
709! Compute nondimensional shape function for diffusion of salinity
710! "Gs1" and its vertical derivative "dGs1dS" evaluated at "hbbl".
711!
712 k_bl=cff_dn*akt(i,j,k,isalt)+cff_up*akt(i,j,k-1,isalt)
713 dk_bl=-cff*(akt(i,j,k,isalt)-akt(i,j,k-1,isalt))
714 gs1(i,j)=k_bl/(zbl*ws(i,j)+eps)
715# ifdef MASKING
716 gs1(i,j)=gs1(i,j)*rmask(i,j)
717# endif
718 dgs1ds(i,j)=min(0.0_r8, k_bl*f1(i,j)-dk_bl/(ws(i,j)+eps))
719# endif
720 END DO
721 END DO
722!
723!-----------------------------------------------------------------------
724! Compute bottom boundary layer mixing coefficients.
725!-----------------------------------------------------------------------
726!
727 DO k=1,n(ng)-1
728 DO j=jstr,jend
729 DO i=istr,iend
730 IF (z_w(i,j,k).lt.hbbl(i,j)) THEN
731!
732! Compute turbulent velocity scales at vertical W-points.
733!
734 depth=z_w(i,j,k)-z_w(i,j,0)
735 IF (bflux(i,j,k).lt.0.0_r8) THEN
736 sigma=min(bl_dpth(i,j),depth)
737 ELSE
738 sigma=depth
739 END IF
740 ustar3=ustar(i,j)*ustar(i,j)*ustar(i,j)
741 zetahat=vonkar*sigma*bflux(i,j,k)
742 zetapar=zetahat/(ustar3+small)
743 IF (zetahat.ge.0.0_r8) THEN ! stable
744 wm(i,j)=vonkar*ustar(i,j)/(1.0_r8+5.0_r8*zetapar)
745 ws(i,j)=wm(i,j)
746 ELSE ! unstable
747 IF (zetapar.gt.lmd_zetam) THEN
748 wm(i,j)=vonkar*ustar(i,j)* &
749 & (1.0_r8-16.0_r8*zetapar)**0.25_r8
750 ELSE
751 wm(i,j)=vonkar*(lmd_am*ustar3-lmd_cm*zetahat)**r3
752 END IF
753 IF (zetapar.gt.lmd_zetas) THEN
754 ws(i,j)=vonkar*ustar(i,j)* &
755 & (1.0_r8-16.0_r8*zetapar)**0.5_r8
756 ELSE
757 ws(i,j)=vonkar*(lmd_as*ustar3-lmd_cs*zetahat)**r3
758 END IF
759 END IF
760!
761! Set polynomial coefficients for shape function.
762!
763 sigma=depth/(hbbl(i,j)-z_w(i,j,0)+eps)
764# ifdef MASKING
765 sigma=sigma*rmask(i,j)
766# endif
767 a1=sigma-2.0_r8
768 a2=3.0_r8-2.0_r8*sigma
769 a3=sigma-1.0_r8
770!
771! Compute nondimesional shape functions.
772!
773 gm=a1+a2*gm1(i,j)+a3*dgm1ds(i,j)
774 gt=a1+a2*gt1(i,j)+a3*dgt1ds(i,j)
775# ifdef SALINITY
776 gs=a1+a2*gs1(i,j)+a3*dgs1ds(i,j)
777# endif
778!
779! Compute boundary layer mixing coefficients, combine them
780! with interior mixing coefficients. Take the maximum estimate
781! of vertical mixing only where the surface and bottom boundary
782! layers overlap. (Do not let the interior overwrite the boundary
783! layer estimate).
784!
785 IF (k.gt.ksbl(i,j)) THEN
786 akv(i,j,k)=max(akv(i,j,k), &
787 & depth*wm(i,j)*(1.0_r8+sigma*gm))
788 akt(i,j,k,itemp)=max(akt(i,j,k,itemp), &
789 & depth*ws(i,j)*(1.0_r8+sigma*gt))
790# ifdef SALINITY
791 akt(i,j,k,isalt)=max(akt(i,j,k,isalt), &
792 & depth*ws(i,j)*(1.0_r8+sigma*gs))
793# endif
794 ELSE
795 akv(i,j,k)=depth*wm(i,j)*(1.0_r8+sigma*gm)
796 akt(i,j,k,itemp)=depth*ws(i,j)*(1.0_r8+sigma*gt)
797# ifdef SALINITY
798 akt(i,j,k,isalt)=depth*ws(i,j)*(1.0_r8+sigma*gs)
799# endif
800 END IF
801 END IF
802 END DO
803 END DO
804 END DO
805
806 RETURN
807 END SUBROUTINE lmd_bkpp_tile
808#endif
809 END MODULE lmd_bkpp_mod
lmd_bkpp_tile
subroutine lmd_bkpp_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, rmask, f, h, hz, z_r, z_w, u, v, pden, srflx, btflx, bustr, bvstr, alpha, beta, bvf, ksbl, akt, akv, kbbl, hbbl)
Definition lmd_bkpp.F:111
lmd_bkpp
subroutine lmd_bkpp(ng, tile)
Definition lmd_bkpp.F:45
lmd_swfrac_tile
subroutine lmd_swfrac_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, zscale, z, swdk)
Definition lmd_swfrac.F:10
bc_2d_mod
Definition bc_2d.F:2
bc_2d_mod::bc_r2d_tile
subroutine bc_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition bc_2d.F:44
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_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_scalars
Definition mod_scalars.F:2
mod_scalars::lmd_betat
real(r8) lmd_betat
Definition mod_scalars.F:1705
mod_scalars::lmd_epsilon
real(r8) lmd_epsilon
Definition mod_scalars.F:1710
mod_scalars::vonkar
real(dp) vonkar
Definition mod_scalars.F:469
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::lmd_cm
real(r8) lmd_cm
Definition mod_scalars.F:1708
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::lmd_cv
real(r8) lmd_cv
Definition mod_scalars.F:1695
mod_scalars::lmd_zetam
real(r8) lmd_zetam
Definition mod_scalars.F:1711
mod_scalars::gorho0
real(dp) gorho0
Definition mod_scalars.F:468
mod_scalars::lmd_zetas
real(r8) lmd_zetas
Definition mod_scalars.F:1712
mod_scalars::isalt
integer isalt
Definition mod_scalars.F:127
mod_scalars::itemp
integer itemp
Definition mod_scalars.F:126
mod_scalars::lmd_ric
real(r8) lmd_ric
Definition mod_scalars.F:1700
mod_scalars::lmd_as
real(r8) lmd_as
Definition mod_scalars.F:1704
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::lmd_am
real(r8) lmd_am
Definition mod_scalars.F:1703
mod_scalars::lmd_cs
real(r8) lmd_cs
Definition mod_scalars.F:1709
mod_scalars::lmd_cekman
real(r8) lmd_cekman
Definition mod_scalars.F:1706
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294
shapiro_mod
Definition shapiro.F:2
shapiro_mod::shapiro2d_tile
subroutine shapiro2d_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, amask, a)
Definition shapiro.F:33