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