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/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpRK3.F90 @ 13328

Last change on this file since 13328 was 13328, checked in by gm, 4 years ago

fixes from last merge qco r12983

File size: 19.7 KB
Line 
1MODULE stpRK3
2   !!======================================================================
3   !!                       ***  MODULE step  ***
4   !! Time-stepping   : manager of the shallow water equation time stepping
5   !!======================================================================
6   !! History :  NEMO !  2020-03  (A. Nasser, G. Madec)  Original code from  4.0.2
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   stpRK3             : Shallow Water time-stepping
11   !!----------------------------------------------------------------------
12   USE step_oce         ! time stepping definition modules
13   USE phycst           ! physical constants
14   USE usrdef_nam
15   !
16   USE iom              ! xIOs server
17   USE domqco
18
19   IMPLICIT NONE
20   PRIVATE
21
22   PUBLIC   stp_RK3   ! called by nemogcm.F90
23   
24   !!----------------------------------------------------------------------
25   !! time level indices
26   !!----------------------------------------------------------------------
27   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !! used by nemo_init
28     
29   !! * Substitutions
30#  include "do_loop_substitute.h90"
31#  include "domzgr_substitute.h90"
32   !!----------------------------------------------------------------------
33   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
34   !! $Id: step.F90 12614 2020-03-26 14:59:52Z gm $
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39#if defined key_agrif
40   RECURSIVE SUBROUTINE stp_RK3( )
41      INTEGER             ::   kstp   ! ocean time-step index
42#else
43   SUBROUTINE stp_RK3( kstp )
44      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
45#endif
46      !!----------------------------------------------------------------------
47      !!                     ***  ROUTINE stp_RK3  ***
48      !!
49      !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)
50      !!
51      !! ** Method  : -1- Update forcings
52      !!              -2- Update the ssh at Naa
53      !!              -3- Compute the momentum trends (Nrhs)
54      !!              -4- Update the horizontal velocity
55      !!              -5- Apply Asselin time filter to uu,vv,ssh
56      !!              -6- Outputs and diagnostics
57      !!----------------------------------------------------------------------
58      INTEGER ::   ji, jj, jk   ! dummy loop indice
59      INTEGER ::   indic        ! error indicator if < 0
60!!gm kcall can be removed, I guess
61      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
62      REAL(wp)::   z1_2rho0,  z5_6,  z3_4  ! local scalars
63     
64      REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars
65      REAL(wp) ::   zve3a, zve3n, zve3b    !   -      -
66      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva
67      !! ---------------------------------------------------------------------
68#if defined key_agrif
69      kstp = nit000 + Agrif_Nb_Step()
70      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
71      IF( lk_agrif_debug ) THEN
72         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
73         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
74      ENDIF
75      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
76# if defined key_iomput
77      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
78# endif
79#endif
80      !
81      IF( ln_timing )   CALL timing_start('stp_RK3')
82      !
83      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
84      ! model timestep
85      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
86      !
87      IF ( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero
88
89      !
90      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
91      ! update I/O and calendar
92      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
93                             indic = 0                ! reset to no error condition
94                             
95      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
96                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
97         IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
98                             CALL iom_init_closedef
99         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
100      ENDIF
101      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
102                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
103      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
104
105      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
106      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
107      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
108      IF( ln_tide    )   CALL tide_update( kstp )                     ! update tide potential
109      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
110      IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries
111                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition
112
113      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
114      ! Ocean physics update
115      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
116      !  LATERAL  PHYSICS
117      !                                                                        ! eddy diffusivity coeff.
118      IF( l_ldfdyn_time )   CALL ldf_dyn( kstp, Nbb )                          ! eddy viscosity coeff.
119
120
121      !======================================================================
122      !======================================================================
123      !                     =====       RK3       =====
124      !======================================================================
125      !======================================================================
126
127     
128      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
129      !  RK3 1st stage Ocean dynamics : hdiv, ssh, e3, u, v, w
130      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
131      rDt   = rn_Dt / 3._wp 
132      r1_Dt = 1._wp / rDt
133     
134                            CALL ssh_nxt       ( kstp, Nbb, Nbb, ssh, Naa )    ! after ssh (includes call to div_hor)
135
136                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
137                         vv(:,:,:,Nrhs) = 0._wp
138
139!!                            CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF) ==> RHS
140 
141                            CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS
142#if defined key_RK3all 
143                            CALL dyn_ldf( kstp, Nbb, Nbb      , uu, vv, Nrhs )  ! lateral mixing
144#endif
145      !
146!!an - calcul du gradient de pression horizontal (explicit)
147      DO_3D_00_00( 1, jpkm1 )
148         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj)
149         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
150      END_3D
151      !
152#if defined key_RK3all 
153      ! add wind stress forcing and layer linear friction to the RHS
154      z5_6 = 5._wp/6._wp
155      DO_3D_00_00(1,jpkm1)
156         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)   &
157            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
158         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)   &
159            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
160      END_3D
161#endif
162!!an
163                            CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit
164      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
165         DO_3D_00_00(1,jpkm1)
166            uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
167            vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
168         END_3D         
169      ELSE
170         DO_3D_00_00(1,jpkm1)       ! flux form : applied on thickness weighted velocity
171            uu(ji,jj,jk,Naa) = (         uu(ji,jj,jk,Nbb )*e3u(ji,jj,jk,Nbb)                              &
172               &                 + rDt * uu(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nbb) * umask(ji,jj,jk)        )   &
173               &               /                           e3t(ji,jj,jk,Naa)
174            vv(ji,jj,jk,Naa) = (         vv(ji,jj,jk,Nbb )*e3v(ji,jj,jk,Nbb)                              &
175               &                 + rDt * vv(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nbb) * vmask(ji,jj,jk)        )   &
176               &               /                           e3t(ji,jj,jk,Naa)
177         END_3D
178      ENDIF
179     
180      CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
181             &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
182
183
184      ! Swap time levels
185      Nrhs= Nnn
186      Nnn = Naa
187      Naa = Nrhs
188
189     
190      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
191      !  RK3 2nd stage Ocean dynamics : hdiv, ssh, e3, u, v, w
192      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
193      rDt   = rn_Dt / 2._wp 
194      r1_Dt = 1._wp / rDt
195     
196                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
197
198                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
199                         vv(:,:,:,Nrhs) = 0._wp
200!!st TBC for dyn_adv
201!!                            CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF) ==> RHS
202 
203                            CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS
204#if defined key_RK3all 
205                            CALL dyn_ldf( kstp, Nbb, Nbb      , uu, vv, Nrhs )  ! lateral mixing
206#endif
207                           
208      !
209!!an - calcul du gradient de pression horizontal (explicit)
210      DO_3D_00_00( 1, jpkm1 )
211         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
212         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
213      END_3D
214      !
215      ! add wind stress forcing and layer linear friction to the RHS
216#if defined key_RK3all
217      z3_4 = 3._wp/4._wp
218      DO_3D_00_00(1,jpkm1)
219         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)   &
220            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
221         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)   &
222            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
223      END_3D
224#endif
225!!an
226                           CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit
227      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
228         DO_3D_00_00(1,jpkm1)
229            uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
230            vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
231         END_3D         
232      ELSE
233         DO_3D_00_00(1,jpkm1)       ! flux form : applied on thickness weighted velocity
234            uu(ji,jj,jk,Naa) = (         uu(ji,jj,jk,Nbb )*e3u(ji,jj,jk,Nbb)                              &
235               &                 + rDt * uu(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nnn) * umask(ji,jj,jk)        )   &
236               &               /                           e3t(ji,jj,jk,Naa)
237            vv(ji,jj,jk,Naa) = (         vv(ji,jj,jk,Nbb )*e3v(ji,jj,jk,Nbb)                              &
238               &                 + rDt * vv(ji,jj,jk,Nrhs)*e3t(ji,jj,jk,Nnn) * vmask(ji,jj,jk)        )   &
239               &               /                           e3t(ji,jj,jk,Naa)
240         END_3D
241      ENDIF
242     
243      CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
244             &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
245
246      ! Swap time levels
247      Nrhs= Nnn
248      Nnn = Naa
249      Naa = Nrhs
250       
251      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
252      !  RK3 3rd stage Ocean dynamics : hdiv, ssh, e3, u, v, w
253      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
254      rDt   = rn_Dt
255      r1_Dt = 1._wp / rDt
256
257                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
258     
259                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
260                         vv(:,:,:,Nrhs) = 0._wp
261
262      IF( ln_bdy     )      CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends
263
264#if defined key_agrif
265      IF(.NOT. Agrif_Root())  & 
266               &            CALL Agrif_Sponge_dyn        ! momentum sponge
267#endif
268!!                            CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF) ==> RHS
269 
270                            CALL dyn_vor( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS
271 
272                            CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
273
274!!an - calcul du gradient de pression horizontal (explicit)
275      DO_3D_00_00( 1, jpkm1 )
276         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
277         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
278      END_3D
279      !
280      ! add wind stress forcing and layer linear friction to the RHS
281      z1_2rho0 = 0.5_wp * r1_rho0
282      DO_3D_00_00(1,jpkm1)
283         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
284            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
285         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
286            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
287      END_3D   
288!!an         
289                            CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit     
290      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
291         DO_3D_11_11(1,jpkm1)
292            zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
293            zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
294            !                                                                  ! Asselin time filter on u,v (Nnn)
295            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)
296            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)
297            !             
298            uu(ji,jj,jk,Naa) = zua
299            vv(ji,jj,jk,Naa) = zva
300         END_3D
301         !
302      ELSE                          ! flux form : applied on thickness weighted velocity
303         DO_3D_11_11(1,jpkm1)
304            zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn)
305            zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn)
306            zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
307            zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
308            !                                                ! LF time stepping
309            zue3a = zue3b + rDt * e3t(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
310            zve3a = zve3b + rDt * e3t(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
311            !
312            uu(ji,jj,jk,Naa) = zue3a / e3t(ji,jj,jk,Naa)   
313            vv(ji,jj,jk,Naa) = zve3a / e3t(ji,jj,jk,Naa)
314         END_3D
315!!st je ne comprends pas l'histoire des e3t et du Nbb et pas du Nnn pour le rhs ?       
316      ENDIF
317
318
319      CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
320             &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
321      ! Swap time levels
322      Nrhs = Nbb
323      Nbb = Naa
324      Naa = Nrhs
325      !
326!                         CALL dom_vvl_sf_update_st( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
327      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
328      ! diagnostics and outputs
329      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
330      IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats
331      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics
332   
333                         CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
334
335      !
336      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
337         
338
339#if defined key_agrif
340      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
341      ! AGRIF
342      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<     
343                         Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
344                         CALL Agrif_Integrate_ChildGrids( stp_RK3 )       ! allows to finish all the Child Grids before updating
345
346                         IF( Agrif_NbStepint() == 0 ) THEN
347                            CALL Agrif_update_all( )                  ! Update all components
348                         ENDIF
349#endif
350      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
351
352      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
353      ! Control
354      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
355                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
356               
357                         
358      IF( kstp == nit000 ) THEN                          ! 1st time step only
359                                        CALL iom_close( numror )   ! close input  ocean restart file
360         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
361         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
362      ENDIF
363
364      !
365#if defined key_iomput
366      IF( kstp == nitend .OR. indic < 0 ) THEN
367                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
368                      IF(lrxios) CALL iom_context_finalize(      crxios_context          )
369      ENDIF
370#endif
371      !
372      IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
373         rDt = 2._wp * rn_Dt   
374         r1_Dt = 1._wp / rDt
375         l_1st_euler = .FALSE.     
376      ENDIF
377      !
378      IF( ln_timing )   CALL timing_stop('stp_RK3')
379      !
380   END SUBROUTINE stp_RK3
381   !
382   !!======================================================================
383END MODULE stpRK3
Note: See TracBrowser for help on using the repository browser.