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