ROMS
Loading...
Searching...
No Matches
gls_prestep.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE gls_prestep_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 predictor step for turbulent kinetic !
13! energy prognostic variables, tke and gls. A NON-conservative, !
14! but constancy preserving, auxiliary advective substep for tke !
15! gls equations is carried out. The result of this substep will !
16! be used to compute advective terms in the corrector substep. !
17! No dissipation terms are included here. !
18! !
19!=======================================================================
20!
21 implicit none
22!
23 PRIVATE
24 PUBLIC :: gls_prestep
25!
26 CONTAINS
27!
28!***********************************************************************
29 SUBROUTINE gls_prestep (ng, tile)
30!***********************************************************************
31!
32 USE mod_param
33 USE mod_grid
34 USE mod_ocean
35 USE mod_mixing
36 USE mod_stepping
37!
38! Imported variable declarations.
39!
40 integer, intent(in) :: ng, tile
41!
42! Local variable declarations.
43!
44 character (len=*), parameter :: myfile = &
45 & __FILE__
46!
47# include "tile.h"
48!
49# ifdef PROFILE
50 CALL wclock_on (ng, inlm, 19, __line__, myfile)
51# endif
52 CALL gls_prestep_tile (ng, tile, &
53 & lbi, ubi, lbj, ubj, &
54 & imins, imaxs, jmins, jmaxs, &
55 & nstp(ng), nnew(ng), &
56# ifdef MASKING
57 & grid(ng) % umask, &
58 & grid(ng) % vmask, &
59# endif
60 & grid(ng) % Huon, &
61 & grid(ng) % Hvom, &
62 & grid(ng) % Hz, &
63 & grid(ng) % pm, &
64 & grid(ng) % pn, &
65 & ocean(ng) % W, &
66# ifdef WEC_VF
67 & ocean(ng) % W_stokes, &
68# endif
69 & mixing(ng) % gls, &
70 & mixing(ng) % tke)
71# ifdef PROFILE
72 CALL wclock_off (ng, inlm, 19, __line__, myfile)
73# endif
74!
75 RETURN
76 END SUBROUTINE gls_prestep
77!
78!***********************************************************************
79 SUBROUTINE gls_prestep_tile (ng, tile, &
80 & LBi, UBi, LBj, UBj, &
81 & IminS, ImaxS, JminS, JmaxS, &
82 & nstp, nnew, &
83# ifdef MASKING
84 & umask, vmask, &
85# endif
86 & Huon, Hvom, Hz, pm, pn, W, &
87# ifdef WEC_VF
88 & W_stokes, &
89# endif
90 & gls, tke)
91!***********************************************************************
92!
93 USE mod_param
94 USE mod_scalars
95!
96 USE exchange_3d_mod, ONLY : exchange_w3d_tile
97# ifdef DISTRIBUTE
98 USE mp_exchange_mod, ONLY : mp_exchange3d
99# endif
100 USE tkebc_mod, ONLY : tkebc_tile
101!
102! Imported variable declarations.
103!
104 integer, intent(in) :: ng, tile
105 integer, intent(in) :: LBi, UBi, LBj, UBj
106 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
107 integer, intent(in) :: nstp, nnew
108!
109# ifdef ASSUMED_SHAPE
110# ifdef MASKING
111 real(r8), intent(in) :: umask(LBi:,LBj:)
112 real(r8), intent(in) :: vmask(LBi:,LBj:)
113# endif
114 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
115 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
116 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
117 real(r8), intent(in) :: pm(LBi:,LBj:)
118 real(r8), intent(in) :: pn(LBi:,LBj:)
119 real(r8), intent(in) :: W(LBi:,LBj:,0:)
120# ifdef WEC_VF
121 real(r8), intent(in) :: W_stokes(LBi:,LBj:,0:)
122# endif
123 real(r8), intent(inout) :: gls(LBi:,LBj:,0:,:)
124 real(r8), intent(inout) :: tke(LBi:,LBj:,0:,:)
125# else
126# ifdef MASKING
127 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
128 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
129# endif
130 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
131 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
132 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
133 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
134 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
135 real(r8), intent(in) :: W(LBi:UBi,LBj:UBj,0:N(ng))
136# ifdef WEC_VF
137 real(r8), intent(in) :: W_stokes(LBi:UBi,LBj:UBj,0:N(ng))
138# endif
139 real(r8), intent(inout) :: gls(LBi:UBi,LBj:UBj,0:N(ng),3)
140 real(r8), intent(inout) :: tke(LBi:UBi,LBj:UBj,0:N(ng),3)
141# endif
142!
143! Local variable declarations.
144!
145 integer :: i, indx, j, k
146
147 real(r8), parameter :: Gamma = 1.0_r8/6.0_r8
148
149 real(r8) :: cff, cff1, cff2, cff3, cff4
150
151 real(r8), dimension(IminS:ImaxS,N(ng)) :: CF
152 real(r8), dimension(IminS:ImaxS,N(ng)) :: FC
153 real(r8), dimension(IminS:ImaxS,N(ng)) :: FCL
154
155 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,N(ng)) :: Hz_half
156
157 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: EF
158 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
159 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FEL
160 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
161 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FXL
162 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: XF
163 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
164 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradL
165
166# include "set_bounds.h"
167!
168!-----------------------------------------------------------------------
169! Predictor step for advection of turbulent kinetic energy variables.
170!-----------------------------------------------------------------------
171!
172! Start computation of auxiliary time step fields tke(:,:,:,n+1/2) and
173! gls(:,:,:,n+1/2) with computation of horizontal advection terms and
174! auxiliary grid-box height field Hz_new()=Hz(:,:,k+1/2,n+1/2);
175! This is effectivey an LF step with subsequent interpolation of the
176! result half step back, using AM3 weights. The LF step and
177! interpolation are perfomed as a single operation, which results in
178! weights cff1,cff2,cff3 below.
179!
180! Either centered fourth-order accurate or standard second order
181! accurate versions are supported.
182!
183! At the same time prepare for corrector step for tke,gls: set tke,
184! gls(:,:,:,nnew) to tke,gls(:,:,:,nstp) multiplied by the
185! corresponding grid-box height. This needs done at this time because
186! array Hz(:,:,:) will overwritten after 2D time stepping with the
187! values computed from zeta(:,:,n+1) rather than zeta(:,:,n), so that
188! the old-time-step Hz will be no longer awailable.
189!
190
191 DO k=1,n(ng)-1
192# ifdef K_C2ADVECTION
193!
194! Second-order, centered differences advection.
195!
196 DO j=jstr,jend
197 DO i=istr,iend+1
198 xf(i,j)=0.5_r8*(huon(i,j,k)+huon(i,j,k+1))
199 fx(i,j)=xf(i,j)* &
200 & 0.5_r8*(tke(i,j,k,nstp)+tke(i-1,j,k,nstp))
201 fxl(i,j)=xf(i,j)* &
202 & 0.5_r8*(gls(i,j,k,nstp)+gls(i-1,j,k,nstp))
203 END DO
204 END DO
205 DO j=jstr,jend+1
206 DO i=istr,iend
207 ef(i,j)=0.5*(hvom(i,j,k)+hvom(i,j,k+1))
208 fe(i,j)=ef(i,j)* &
209 & 0.5*(tke(i,j,k,nstp)+tke(i,j-1,k,nstp))
210 fel(i,j)=ef(i,j)* &
211 & 0.5*(gls(i,j,k,nstp)+gls(i,j-1,k,nstp))
212 END DO
213 END DO
214# else
215!
216! Fourth-order, centered differences advection.
217!
218 DO j=jstr,jend
219 DO i=istrm1,iendp2
220 grad(i,j)=(tke(i,j,k,nstp)-tke(i-1,j,k,nstp))
221# ifdef MASKING
222 grad(i,j)=grad(i,j)*umask(i,j)
223# endif
224 gradl(i,j)=(gls(i,j,k,nstp)-gls(i-1,j,k,nstp))
225# ifdef MASKING
226 gradl(i,j)=gradl(i,j)*umask(i,j)
227# endif
228 END DO
229 END DO
230 IF (.not.(compositegrid(iwest,ng).or.ewperiodic(ng))) THEN
231 IF (domain(ng)%Western_Edge(tile)) THEN
232 DO j=jstr,jend
233 grad(istr-1,j)=grad(istr,j)
234 gradl(istr-1,j)=gradl(istr,j)
235 END DO
236 END IF
237 END IF
238 IF (.not.(compositegrid(ieast,ng).or.ewperiodic(ng))) THEN
239 IF (domain(ng)%Eastern_Edge(tile)) THEN
240 DO j=jstr,jend
241 grad(iend+2,j)=grad(iend+1,j)
242 gradl(iend+2,j)=gradl(iend+1,j)
243 END DO
244 END IF
245 END IF
246 cff=1.0_r8/6.0_r8
247 DO j=jstr,jend
248 DO i=istr,iend+1
249 xf(i,j)=0.5_r8*(huon(i,j,k)+huon(i,j,k+1))
250 fx(i,j)=xf(i,j)* &
251 & 0.5_r8*(tke(i-1,j,k,nstp)+tke(i,j,k,nstp)- &
252 & cff*(grad(i+1,j)-grad(i-1,j)))
253 fxl(i,j)=xf(i,j)* &
254 & 0.5_r8*(gls(i-1,j,k,nstp)+gls(i,j,k,nstp)- &
255 & cff*(gradl(i+1,j)-gradl(i-1,j)))
256 END DO
257 END DO
258!
259 DO j=jstrm1,jendp2
260 DO i=istr,iend
261 grad(i,j)=(tke(i,j,k,nstp)-tke(i,j-1,k,nstp))
262# ifdef MASKING
263 grad(i,j)=grad(i,j)*vmask(i,j)
264# endif
265 gradl(i,j)=(gls(i,j,k,nstp)-gls(i,j-1,k,nstp))
266# ifdef MASKING
267 gradl(i,j)=gradl(i,j)*vmask(i,j)
268# endif
269 END DO
270 END DO
271 IF (.not.(compositegrid(isouth,ng).or.nsperiodic(ng))) THEN
272 IF (domain(ng)%Southern_Edge(tile)) THEN
273 DO i=istr,iend
274 grad(i,jstr-1)=grad(i,jstr)
275 gradl(i,jstr-1)=gradl(i,jstr)
276 END DO
277 END IF
278 END IF
279 IF (.not.(compositegrid(inorth,ng).or.nsperiodic(ng))) THEN
280 IF (domain(ng)%Northern_Edge(tile)) THEN
281 DO i=istr,iend
282 grad(i,jend+2)=grad(i,jend+1)
283 gradl(i,jend+2)=gradl(i,jend+1)
284 END DO
285 END IF
286 END IF
287 cff=1.0_r8/6.0_r8
288 DO j=jstr,jend+1
289 DO i=istr,iend
290 ef(i,j)=0.5_r8*(hvom(i,j,k)+hvom(i,j,k+1))
291 fe(i,j)=ef(i,j)* &
292 & 0.5_r8*(tke(i,j-1,k,nstp)+tke(i,j,k,nstp)- &
293 & cff*(grad(i,j+1)-grad(i,j-1)))
294 fel(i,j)=ef(i,j)* &
295 & 0.5_r8*(gls(i,j-1,k,nstp)+gls(i,j,k,nstp)- &
296 & cff*(gradl(i,j+1)-gradl(i,j-1)))
297 END DO
298 END DO
299# endif
300!
301! Time-step horizontal advection.
302!
303 IF (iic(ng).eq.ntfirst(ng)) THEN
304 cff1=1.0_r8
305 cff2=0.0_r8
306 cff3=0.5_r8*dt(ng)
307 indx=nstp
308 ELSE
309 cff1=0.5_r8+gamma
310 cff2=0.5_r8-gamma
311 cff3=(1.0_r8-gamma)*dt(ng)
312 indx=3-nstp
313 END IF
314 DO j=jstr,jend
315 DO i=istr,iend
316 cff=0.5_r8*(hz(i,j,k)+hz(i,j,k+1))
317 cff4=cff3*pm(i,j)*pn(i,j)
318 hz_half(i,j,k)=cff-cff4*(xf(i+1,j)-xf(i,j)+ &
319 & ef(i,j+1)-ef(i,j))
320 tke(i,j,k,3)=cff*(cff1*tke(i,j,k,nstp)+ &
321 & cff2*tke(i,j,k,indx))- &
322 & cff4*(fx(i+1,j)-fx(i,j)+ &
323 & fe(i,j+1)-fe(i,j))
324 gls(i,j,k,3)=cff*(cff1*gls(i,j,k,nstp)+ &
325 & cff2*gls(i,j,k,indx))- &
326 & cff4*(fxl(i+1,j)-fxl(i,j)+ &
327 & fel(i,j+1)-fel(i,j))
328 tke(i,j,k,nnew)=cff*tke(i,j,k,nstp)
329 gls(i,j,k,nnew)=cff*gls(i,j,k,nstp)
330 END DO
331 END DO
332 END DO
333!
334! Compute vertical advection term.
335!
336 DO j=jstr,jend
337# ifdef K_C2ADVECTION
338 DO k=1,n(ng)
339 DO i=istr,iend
340 cf(i,k)=0.5_r8*(w(i,j,k)+w(i,j,k-1))
341# ifdef WEC_VF
342 cf(i,k)=cf(i,k)+0.5_r8*(w_stokes(i,j,k)+w_stokes(i,j,k-1))
343# endif
344 fc(i,k)=cf(i,k)* &
345 & 0.5_r8*(tke(i,j,k-1,nstp)+tke(i,j,k,nstp))
346 fcl(i,k)=cf(i,k)* &
347 & 0.5_r8*(gls(i,j,k-1,nstp)+gls(i,j,k,nstp))
348 END DO
349 END DO
350# else
351 cff1=7.0_r8/12.0_r8
352 cff2=1.0_r8/12.0_r8
353 DO k=2,n(ng)-1
354 DO i=istr,iend
355 cf(i,k)=0.5_r8*(w(i,j,k)+w(i,j,k-1))
356# ifdef WEC_VF
357 cf(i,k)=cf(i,k)+0.5_r8*(w_stokes(i,j,k)+w_stokes(i,j,k-1))
358# endif
359 fc(i,k)=cf(i,k)*(cff1*(tke(i,j,k-1,nstp)+ &
360 & tke(i,j,k ,nstp))- &
361 & cff2*(tke(i,j,k-2,nstp)+ &
362 & tke(i,j,k+1,nstp)))
363 fcl(i,k)=cf(i,k)*(cff1*(gls(i,j,k-1,nstp)+ &
364 & gls(i,j,k ,nstp))- &
365 & cff2*(gls(i,j,k-2,nstp)+ &
366 & gls(i,j,k+1,nstp)))
367 END DO
368 END DO
369 cff1=1.0_r8/3.0_r8
370 cff2=5.0_r8/6.0_r8
371 cff3=1.0_r8/6.0_r8
372 DO i=istr,iend
373 cf(i,1)=0.5*(w(i,j,0)+w(i,j,1))
374# ifdef WEC_VF
375 cf(i,1)=cf(i,1)+0.5_r8*(w_stokes(i,j,0)+w_stokes(i,j,1))
376# endif
377 fc(i,1)=cf(i,1)*(cff1*tke(i,j,0,nstp)+ &
378 & cff2*tke(i,j,1,nstp)- &
379 & cff3*tke(i,j,2,nstp))
380 fcl(i,1)=cf(i,1)*(cff1*gls(i,j,0,nstp)+ &
381 & cff2*gls(i,j,1,nstp)- &
382 & cff3*gls(i,j,2,nstp))
383 cf(i,n(ng))=0.5*(w(i,j,n(ng))+w(i,j,n(ng)-1))
384# ifdef WEC_VF
385 cf(i,n(ng))=cf(i,n(ng))+0.5_r8* &
386 & (w_stokes(i,j,n(ng))+w_stokes(i,j,n(ng)-1))
387# endif
388 fc(i,n(ng))=cf(i,n(ng))*(cff1*tke(i,j,n(ng) ,nstp)+ &
389 & cff2*tke(i,j,n(ng)-1,nstp)- &
390 & cff3*tke(i,j,n(ng)-2,nstp))
391 fcl(i,n(ng))=cf(i,n(ng))*(cff1*gls(i,j,n(ng) ,nstp)+ &
392 & cff2*gls(i,j,n(ng)-1,nstp)- &
393 & cff3*gls(i,j,n(ng)-2,nstp))
394 END DO
395# endif
396!
397! Time-step vertical advection term.
398!
399 IF (iic(ng).eq.ntfirst(ng)) THEN
400 cff3=0.5_r8*dt(ng)
401 ELSE
402 cff3=(1.0_r8-gamma)*dt(ng)
403 END IF
404 DO k=1,n(ng)-1
405 DO i=istr,iend
406 cff4=cff3*pm(i,j)*pn(i,j)
407 hz_half(i,j,k)=hz_half(i,j,k)-cff4*(cf(i,k+1)-cf(i,k))
408 cff1=1.0_r8/hz_half(i,j,k)
409 tke(i,j,k,3)=cff1*(tke(i,j,k,3)- &
410 & cff4*(fc(i,k+1)-fc(i,k)))
411 gls(i,j,k,3)=cff1*(gls(i,j,k,3)- &
412 & cff4*(fcl(i,k+1)-fcl(i,k)))
413 END DO
414 END DO
415 END DO
416!
417! Apply lateral boundary conditions.
418!
419 CALL tkebc_tile (ng, tile, &
420 & lbi, ubi, lbj, ubj, n(ng), &
421 & imins, imaxs, jmins, jmaxs, &
422 & 3, nstp, &
423 & gls, tke)
424
425 IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
426 CALL exchange_w3d_tile (ng, tile, &
427 & lbi, ubi, lbj, ubj, 0, n(ng), &
428 & tke(:,:,:,3))
429 CALL exchange_w3d_tile (ng, tile, &
430 & lbi, ubi, lbj, ubj, 0, n(ng), &
431 & gls(:,:,:,3))
432 END IF
433
434# ifdef DISTRIBUTE
435 CALL mp_exchange3d (ng, tile, inlm, 2, &
436 & lbi, ubi, lbj, ubj, 0, n(ng), &
437 & nghostpoints, &
438 & ewperiodic(ng), nsperiodic(ng), &
439 & tke(:,:,:,3), &
440 & gls(:,:,:,3))
441# endif
442!
443 RETURN
444 END SUBROUTINE gls_prestep_tile
445#endif
446 END MODULE gls_prestep_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_prestep_mod
Definition gls_prestep.F:2
gls_prestep_mod::gls_prestep_tile
subroutine gls_prestep_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, nnew, umask, vmask, huon, hvom, hz, pm, pn, w, gls, tke)
Definition gls_prestep.F:91
gls_prestep_mod::gls_prestep
subroutine, public gls_prestep(ng, tile)
Definition gls_prestep.F:30
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_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_scalars
Definition mod_scalars.F:2
mod_scalars::iic
integer, dimension(:), allocatable iic
Definition mod_scalars.F:249
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
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::compositegrid
logical, dimension(:,:), allocatable compositegrid
Definition mod_scalars.F:493
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::ntfirst
integer, dimension(:), allocatable ntfirst
Definition mod_scalars.F:377
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_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