Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
MODULE stprk3
!!======================================================================
!! *** MODULE stprk3 ***
!! Time-stepping : manager of the shallow water equation time stepping
!! 3rd order Runge-Kutta scheme
!!======================================================================
!! History : NEMO ! 2020-03 (A. Nasser, G. Madec) Original code from 4.0.2
!! - ! 2020-10 (S. Techene, G. Madec) cleanning
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! stp_RK3 : RK3 Shallow Water Eq. time-stepping
!!----------------------------------------------------------------------
USE stp_oce ! modules used in nemo_init and stp_RK3
!
USE domqco ! quasi-eulerian coordinate
USE phycst ! physical constants
USE usrdef_nam ! user defined namelist parameters
IMPLICIT NONE
PRIVATE
PUBLIC stp_RK3 ! called by nemogcm.F90
! !** time level indices **!
INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs !: used by nemo_init
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: step.F90 12614 2020-03-26 14:59:52Z gm $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE stp_RK3( kstp )
!!----------------------------------------------------------------------
!! *** ROUTINE stp_RK3 ***
!!
!! ** Purpose : - RK3 Time stepping scheme for shallow water Eq.
!!
!! ** Method : 3rd order time stepping scheme which has 3 stages
!! * Update calendar and forcings
!! stage 1 : n ==> n+1/3 using variables at n
!! - Compute the rhs of momentum
!! - Time step ssh at Naa (n+1/3)
!! - Time step u,v at Naa (n+1/3)
!! - Swap time indices
!! stage 2 : n ==> n+1/2 using variables at n and n+1/3
!! - Compute the rhs of momentum
!! - Time step ssh at Naa (n+1/2)
!! - Time step u,v at Naa (n+1/2)
!! - Swap time indices
!! stage 3 : n ==> n+1 using variables at n and n+1/2
!! - Compute the rhs of momentum
!! - Time step ssh at Naa (n+1)
!! - Time step u,v at Naa (n+1)
!! - Swap time indices
!! * Outputs and diagnostics
!!
!! NB: in stages 1 and 2 lateral mixing and forcing are not taken
!! into account in the momentum RHS execpt if key_RK3all is used
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: kstp ! ocean time-step index
!
INTEGER :: ji, jj, jk ! dummy loop indice
INTEGER :: indic ! error indicator if < 0
REAL(wp):: z1_2rho0, z5_6, z3_4 ! local scalars
REAL(wp):: zue3a, zue3b, zua, zrhs_u ! local scalars
REAL(wp):: zve3a, zve3b, zva, zrhs_v ! - -
!! ---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('stp_RK3')
!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! model timestep
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
IF ( kstp == nit000 ) ww(:,:,:) = 0._wp ! initialize vertical velocity one for all to zero
!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! update I/O and calendar
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
indic = 0 ! reset to no error condition
IF( kstp == nit000 ) THEN ! initialize IOM context
CALL iom_init( cxios_context, ld_closedef=.FALSE. ) ! for model grid (including possible AGRIF zoom)
CALL iom_init_closedef
ENDIF
IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init)
CALL iom_setkt( kstp - nit000 + 1, cxios_context ) ! tell IOM we are at time step kstp
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Update external forcing (SWE: surface boundary condition only)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Ocean physics update (SWE: eddy viscosity only)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff.
!======================================================================
!======================================================================
! ===== RK3 =====
!======================================================================
!======================================================================
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RK3 1st stage Ocean dynamics : u, v, ssh
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
rDt = rn_Dt / 3._wp
r1_Dt = 1._wp / rDt
!
! !== RHS of the momentum Eq. ==!
!
uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero
vv(:,:,:,Nrhs) = 0._wp
CALL dyn_adv( kstp, Nbb, Nbb, uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS
CALL dyn_vor( kstp, Nbb, uu, vv, Nrhs ) ! vorticity ==> RHS
#if defined key_RK3all
CALL dyn_ldf( kstp, Nbb, Nbb, uu, vv, Nrhs ) ! lateral mixing
#endif
z5_6 = 5._wp/6._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
! ! horizontal pressure gradient
zrhs_u = - grav * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj)
zrhs_v = - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
#if defined key_RK3all
! ! wind stress and layer friction
zrhs_u = zrhs_u + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utauU(ji,jj) ) / e3u(ji,jj,jk,Nbb) &
zrhs_v = zrhs_v + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtauV(ji,jj) ) / e3v(ji,jj,jk,Nbb) &
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
& - rn_rfr * vv(ji,jj,jk,Nbb)
#endif
! ! ==> RHS
uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
END_3D
!
! !== Time stepping of ssh Eq. ==! (and update r3_Naa)
!
CALL ssh_nxt( kstp, Nbb, Nbb, ssh, Naa ) ! after ssh
! ! after ssh/h_0 ratio
CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
!
! !== Time stepping of momentum Eq. ==!
!
IF( ln_dynadv_vec ) THEN ! vector invariant form : applied on velocity
DO_3D( 0, 0, 0, 0, 1,jpkm1)
uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
END_3D
ELSE
DO_3D( 0, 0, 0, 0, 1,jpkm1) ! flux form : applied on thickness weighted velocity
zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
zue3a = zue3b + rDt * e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
zve3a = zve3b + rDt * e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
!
uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)
vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
END_3D
ENDIF
!
CALL lbc_lnk( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
IF (nn_hls==2) CALL lbc_lnk( 'stp_RK3', r3u(:,:,Naa), 'U', 1., r3v(:,:,Naa), 'V', 1.)
!
! !== Swap time levels ==!
Nrhs= Nnn
Nnn = Naa
Naa = Nrhs
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RK3 2nd stage Ocean dynamics : hdiv, ssh, e3, u, v, w
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
rDt = rn_Dt / 2._wp
r1_Dt = 1._wp / rDt
!
! !== RHS of the momentum Eq. ==!
!
uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero
vv(:,:,:,Nrhs) = 0._wp
CALL dyn_adv( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS
CALL dyn_vor( kstp, Nnn, uu, vv, Nrhs ) ! vorticity ==> RHS
#if defined key_RK3all
CALL dyn_ldf( kstp, Nbb, Nbb, uu, vv, Nrhs ) ! lateral mixing
#endif
!
z3_4 = 3._wp/4._wp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
! ! horizontal pressure gradient
zrhs_u = - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
zrhs_v = - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
#if defined key_RK3all
! ! wind stress and layer friction
zrhs_u = zrhs_u + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utauU(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
zrhs_v = zrhs_v + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtauV(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
& - rn_rfr * vv(ji,jj,jk,Nbb)
#endif
! ! ==> RHS
uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
END_3D
!
! !== Time stepping of ssh Eq. ==! (and update r3_Naa)
!
CALL ssh_nxt( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh
! ! after ssh/h_0 ratio
CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
!
! !== Time stepping of momentum Eq. ==!
!
IF( ln_dynadv_vec ) THEN ! vector invariant form : applied on velocity
DO_3D( 0, 0, 0, 0, 1,jpkm1)
uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
END_3D
ELSE
DO_3D( 0, 0, 0, 0, 1,jpkm1) ! flux form : applied on thickness weighted velocity
zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
zue3a = zue3b + rDt * e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
zve3a = zve3b + rDt * e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
!
uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)
vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
END_3D
ENDIF
!
CALL lbc_lnk( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
IF (nn_hls==2) CALL lbc_lnk( 'stp_RK3', r3u(:,:,Naa), 'U', 1., r3v(:,:,Naa), 'V', 1.)
!
! !== Swap time levels ==!
Nrhs= Nnn
Nnn = Naa
Naa = Nrhs
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RK3 3rd stage Ocean dynamics : hdiv, ssh, e3, u, v, w
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
rDt = rn_Dt
r1_Dt = 1._wp / rDt
!
! !== RHS of the momentum Eq. ==!
!
uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero
vv(:,:,:,Nrhs) = 0._wp
!
CALL dyn_adv( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS
CALL dyn_vor( kstp, Nnn, uu, vv, Nrhs ) ! vorticity ==> RHS
CALL dyn_ldf( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! lateral mixing
z1_2rho0 = 0.5_wp * r1_rho0
DO_3D( 0, 0, 0, 0, 1,jpkm1 )
! ! horizontal pressure gradient
zrhs_u = - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
zrhs_v = - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
! ! wind stress and layer friction
zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utauU(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
& - rn_rfr * vv(ji,jj,jk,Nbb)
! ! ==> RHS
uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
END_3D
!
! !== Time stepping of ssh Eq. ==! (and update r3_Naa)
!
CALL ssh_nxt( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh
! ! after ssh/h_0 ratio
CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
!
! !== Time stepping of momentum Eq. ==!
!
IF( ln_dynadv_vec ) THEN ! vector invariant form : applied on velocity
DO_3D( 0, 0, 0, 0, 1,jpkm1)
uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
END_3D
!
ELSE ! flux form : applied on thickness weighted velocity
DO_3D( 0, 0, 0, 0, 1,jpkm1)
zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
zue3a = zue3b + rDt * e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
zve3a = zve3b + rDt * e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
!
uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)
vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
END_3D
ENDIF
!
CALL lbc_lnk( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
IF (nn_hls==2) CALL lbc_lnk( 'stp_RK3', r3u(:,:,Naa), 'U', 1., r3v(:,:,Naa), 'V', 1.)
!
! !== Swap time levels ==!
!
Nrhs = Nbb
Nbb = Naa
Naa = Nrhs
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! diagnostics and outputs at Nbb (i.e. the just computed time step)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( ln_diacfl ) CALL dia_cfl ( kstp, Nbb ) ! Courant number diagnostics
CALL dia_wri ( kstp, Nbb ) ! ocean model: outputs
!
IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nbb ) ! write output ocean restart file
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Control
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL stp_ctl ( kstp , Nbb )
IF( kstp == nit000 ) THEN ! 1st time step only
CALL iom_close( numror ) ! close input ocean restart file
IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce
IF(lwm .AND. numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice (if exist)
ENDIF
!
#if defined key_xios
IF( kstp == nitend .OR. indic < 0 ) THEN
CALL iom_context_finalize( cxios_context )
ENDIF
#endif
!
IF( ln_timing ) CALL timing_stop('stp_RK3')
!
END SUBROUTINE stp_RK3
!!======================================================================
END MODULE stprk3