New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
stprk3.F90 in NEMO/trunk/src/SWE – NEMO

source: NEMO/trunk/src/SWE/stprk3.F90 @ 14143

Last change on this file since 14143 was 14143, checked in by techene, 4 years ago

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