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
MODULE stpmlf
!!======================================================================
!! *** MODULE stpMLF ***
!! Time-stepping : manager of the shallow water equation time stepping
!! Modified Leap Frog scheme
!!======================================================================
!! History : NEMO ! 2020-03 (A. Nasser, G. Madec) Original code from 4.0.2
!! - ! 2020-10 (S. Techene, G. Madec) cleanning
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! stp_MLF : MLF Shallow Water Eq. time-stepping
!!----------------------------------------------------------------------
USE stp_oce ! modules used in nemo_init and stp_MLF
!
USE domqco ! quasi-eulerian coordinate
USE phycst ! physical constants
USE usrdef_nam ! user defined namelist parameters
!!st USE usrdef_sbc ! user defined surface boundary cond
IMPLICIT NONE
PRIVATE
PUBLIC stp_MLF ! 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_MLF( kstp )
!!----------------------------------------------------------------------
!! *** ROUTINE stp_MLF ***
!!
!! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)
!!
!! ** Method : -1- Update forcings
!! -2- Update the ssh at Naa
!! -3- Compute the momentum trends (Nrhs)
!! -4- Update the horizontal velocity
!! -5- Apply Asselin time filter to uu,vv,ssh
!! -6- Outputs and diagnostics
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kstp ! ocean time-step index
!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: indic ! error indicator if < 0
REAL(wp):: z1_2rho0 ! local scalars
REAL(wp):: zrhs_u, zue3a, zue3n, zue3b, zua ! local scalars
REAL(wp):: zrhs_v, zve3a, zve3n, zve3b, zva ! - -
!! ---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('stp_MLF')
!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! model timestep
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
IF( l_1st_euler ) THEN ! start or restart with Euler 1st time-step
rDt = rn_Dt
r1_Dt = 1._wp / rDt
ENDIF
IF( kstp == nit000 ) ww(:,:,:) = 0._wp ! initialize vertical velocity once for all to zero
!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! update I/O and calendar
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
indic = 0 ! reset to no error condition
IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
CALL iom_init( cxios_context, ld_closedef=.FALSE. ) ! for model grid (including passible 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.
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RHS of horizontal velocity
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
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 ==> RHS
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) &
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
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
204
205
206
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
& - 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)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! ! Leap Frog time stepping ==> ssh_Naa and r3_Naa
CALL ssh_nxt( kstp, Nbb, Nnn, ssh , Naa ) ! after ssh
! ! after ssh/h_0 ratio explicit
CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
! ! Asselin filter ==> ssh_Nnn filtered
IF ( .NOT.( l_1st_euler ) ) THEN ! Time filtering of now ssh
ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2._wp * ssh(:,:,Nnn) + ssh(:,:,Naa) )
ENDIF
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Time stepping of dynamics (u,v)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( ln_dynadv_vec ) THEN ! vector invariant form : applied on velocity
IF( l_1st_euler ) THEN ! Euler time stepping (no Asselin filter)
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 ! Leap Frog time stepping + Asselin filter
DO_3D( 0,0, 0,0, 1,jpkm1)
zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
! ! Asselin time filter on u,v (Nnn)
uu(ji,jj,jk,Nnn) = uu(ji,jj,jk,Nnn) + rn_atfp * (uu(ji,jj,jk,Nbb) - 2._wp * uu(ji,jj,jk,Nnn) + zua)
vv(ji,jj,jk,Nnn) = vv(ji,jj,jk,Nnn) + rn_atfp * (vv(ji,jj,jk,Nbb) - 2._wp * vv(ji,jj,jk,Nnn) + zva)
!
uu(ji,jj,jk,Naa) = zua
vv(ji,jj,jk,Naa) = zva
END_3D
! ! Update r3_Nnn
CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn) ) ! now ssh/h_0 ratio from filtrered ssh
ENDIF
!
ELSE ! flux form : applied on thickness weighted velocity
IF( l_1st_euler ) THEN ! Euler time stepping (no Asselin filter)
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)
! ! Euler time stepping
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
ELSE ! Leap Frog time stepping + Asselin filter
CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:) ) ! now ssh/h_0 ratio from filtrered ssh
DO_3D( 0,0, 0,0, 1,jpkm1)
zue3n = ( 1._wp + r3u(ji,jj,Nnn) ) * uu(ji,jj,jk,Nnn)
zve3n = ( 1._wp + r3v(ji,jj,Nnn) ) * vv(ji,jj,jk,Nnn)
zue3b = ( 1._wp + r3u(ji,jj,Nbb) ) * uu(ji,jj,jk,Nbb)
zve3b = ( 1._wp + r3v(ji,jj,Nbb) ) * vv(ji,jj,jk,Nbb)
! ! LF time stepping
zue3a = zue3b + rDt * ( 1._wp + r3u(ji,jj,Nnn) ) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
zve3a = zve3b + rDt * ( 1._wp + r3v(ji,jj,Nnn) ) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
! ! Asselin time filter on u,v (Nnn)
uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ( 1._wp + r3u_f(ji,jj) )
vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ( 1._wp + r3v_f(ji,jj) )
!
uu(ji,jj,jk,Naa) = zue3a / ( 1._wp + r3u(ji,jj,Naa) )
vv(ji,jj,jk,Naa) = zve3a / ( 1._wp + r3v(ji,jj,Naa) )
END_3D
! ! Update r3_Nnn with time filtered values
r3t(:,:,Nnn) = r3t_f(:,:)
r3u(:,:,Nnn) = r3u_f(:,:)
r3v(:,:,Nnn) = r3v_f(:,:)
ENDIF
ENDIF
CALL lbc_lnk( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1., & !* local domain boundaries
& uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
IF (nn_hls==2) CALL lbc_lnk( 'stp_MLF', r3u(:,:,Naa), 'U', 1., r3v(:,:,Naa), 'V', 1.)
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Set boundary conditions, time filter and swap time levels
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! Swap time levels
Nrhs = Nbb
Nbb = Nnn
Nnn = Naa
Naa = Nrhs
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! diagnostics and outputs
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics
CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs
!
IF( lrst_oce ) CALL rst_write( kstp, Nbb, Nnn ) ! write output ocean restart file
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Control
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
CALL stp_ctl ( kstp, Nnn )
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
!!st : cxios_context needed ? because opened earlier ???
CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
ENDIF
#endif
!
IF( l_1st_euler ) THEN ! recover Leap-frog timestep
rDt = 2._wp * rn_Dt
r1_Dt = 1._wp / rDt
l_1st_euler = .FALSE.
ENDIF
!
IF( ln_timing ) CALL timing_stop('stp_MLF')
!
END SUBROUTINE stp_MLF
!!======================================================================
END MODULE stpmlf