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

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90 @ 12614

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

first Shallow Water Eq. update
  • Property svn:keywords set to Id
File size: 16.6 KB
Line 
1MODULE step
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   !!   stp             : 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
18   IMPLICIT NONE
19   PRIVATE
20
21   PUBLIC   stp   ! called by nemogcm.F90
22   
23   !!----------------------------------------------------------------------
24   !! time level indices
25   !!----------------------------------------------------------------------
26   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !! used by nemo_init
27     
28   !! * Substitutions
29#  include "do_loop_substitute.h90"
30   !!----------------------------------------------------------------------
31   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
32   !! $Id$
33   !! Software governed by the CeCILL license (see ./LICENSE)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37#if defined key_agrif
38   RECURSIVE SUBROUTINE stp( )
39      INTEGER             ::   kstp   ! ocean time-step index
40#else
41   SUBROUTINE stp( kstp )
42      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
43#endif
44      !!----------------------------------------------------------------------
45      !!                     ***  ROUTINE stp  ***
46      !!
47      !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)
48      !!
49      !! ** Method  : -1- Update forcings
50      !!              -2- Update the ssh at Naa
51      !!              -3- Compute the momentum trends (Nrhs)
52      !!              -4- Update the horizontal velocity
53      !!              -5- Apply Asselin time filter to uu,vv,ssh
54      !!              -6- Outputs and diagnostics
55      !!----------------------------------------------------------------------
56      INTEGER ::   ji, jj, jk   ! dummy loop indice
57      INTEGER ::   indic        ! error indicator if < 0
58!!gm kcall can be removed, I guess
59      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
60      REAL(wp)::   z1_2rho0     ! local scalars
61     
62      REAL(wp) ::   zue3a, zue3n, zue3b    ! local scalars
63      REAL(wp) ::   zve3a, zve3n, zve3b    !   -      -
64      REAL(wp) ::   ze3t_tf, ze3u_tf, ze3v_tf, zua, zva
65      !! ---------------------------------------------------------------------
66#if defined key_agrif
67      kstp = nit000 + Agrif_Nb_Step()
68      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
69      IF( lk_agrif_debug ) THEN
70         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
71         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
72      ENDIF
73      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
74# if defined key_iomput
75      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
76# endif
77#endif
78      !
79      IF( ln_timing )   CALL timing_start('stp')
80      !
81      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
82      ! model timestep
83      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
84      !
85      IF( l_1st_euler ) THEN     ! start or restart with Euler 1st time-step
86         rDt   =  rn_Dt   
87         r1_Dt = 1._wp / rDt
88      ENDIF
89
90      IF ( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero
91
92      !
93      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
94      ! update I/O and calendar
95      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
96                             indic = 0                ! reset to no error condition
97                             
98      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
99                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
100         IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
101                             CALL iom_init_closedef
102         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
103      ENDIF
104      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
105                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
106      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
107
108      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
109      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
110      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
111      IF( ln_tide    )   CALL tide_update( kstp )                     ! update tide potential
112      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
113      IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries
114                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition
115
116      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
117      ! Ocean physics update
118      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
119      !  LATERAL  PHYSICS
120      !                                                                        ! eddy diffusivity coeff.
121      IF( l_ldfdyn_time                    )   CALL ldf_dyn( kstp, Nbb )       ! eddy viscosity coeff.
122
123      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
124      !  Ocean dynamics : hdiv, ssh, e3, u, v, w
125      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
126
127                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
128
129      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors
130                                                       
131                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
132                         vv(:,:,:,Nrhs) = 0._wp
133
134      IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends
135
136#if defined key_agrif
137      IF(.NOT. Agrif_Root())  & 
138               &         CALL Agrif_Sponge_dyn        ! momentum sponge
139#endif
140                         CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)   ==> RHS
141 
142                         CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity              ==> RHS
143 
144                         CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
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,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
149         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
150      END_3D
151      !
152      ! add wind stress forcing and layer linear friction to the RHS
153      z1_2rho0 = 0.5_wp * r1_rho0
154      DO_3D_00_00(1,jpkm1)
155         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
156            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
157         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
158            &                                  - rn_rfr * vv(ji,jj,jk,Nbb) 
159      END_3D   
160!!an         
161 
162      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
163      ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3
164      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
165     
166!!an futur module dyn_nxt (a la place de dyn_atf)
167     
168      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
169         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter)
170            DO_3D_00_00(1,jpkm1)
171               uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
172               vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
173             END_3D         
174         ELSE                             ! Leap Frog time stepping + Asselin filter         
175            DO_3D_11_11(1,jpkm1)
176               zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
177               zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
178               !                                                                  ! Asselin time filter on u,v (Nnn)
179               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)
180               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)
181               !
182               ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn)  + e3u(ji,jj,jk,Naa) )
183               ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn)  + e3v(ji,jj,jk,Naa) )
184               ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn)  + e3t(ji,jj,jk,Naa) )
185               !
186               e3u(ji,jj,jk,Nnn) = ze3u_tf
187               e3v(ji,jj,jk,Nnn) = ze3v_tf
188               e3t(ji,jj,jk,Nnn) = ze3t_tf
189               !             
190               uu(ji,jj,jk,Naa) = zua
191               vv(ji,jj,jk,Naa) = zva
192             END_3D   
193         ENDIF
194         !
195      ELSE                          ! flux form : applied on thickness weighted velocity
196         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter)
197            DO_3D_00_00(1,jpkm1)
198               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
199               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
200               !                                                ! LF time stepping
201               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
202               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
203               !                                               
204               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)   
205               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)     
206             END_3D         
207         ELSE                             ! Leap Frog time stepping + Asselin filter         
208            DO_3D_11_11(1,jpkm1)
209               zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn)
210               zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn)
211               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
212               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
213               !                                                ! LF time stepping
214               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
215               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
216               !                                                ! Asselin time filter on e3u/v/t
217               ze3u_tf = e3u(ji,jj,jk,Nnn) + rn_atfp * ( e3u(ji,jj,jk,Nbb) - 2._wp * e3u(ji,jj,jk,Nnn)  + e3u(ji,jj,jk,Naa) ) 
218               ze3v_tf = e3v(ji,jj,jk,Nnn) + rn_atfp * ( e3v(ji,jj,jk,Nbb) - 2._wp * e3v(ji,jj,jk,Nnn)  + e3v(ji,jj,jk,Naa) )
219               ze3t_tf = e3t(ji,jj,jk,Nnn) + rn_atfp * ( e3t(ji,jj,jk,Nbb) - 2._wp * e3t(ji,jj,jk,Nnn)  + e3t(ji,jj,jk,Naa) ) 
220               !                                                ! Asselin time filter on u,v (Nnn)
221               uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_tf
222               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_tf           
223               !
224               e3u(ji,jj,jk,Nnn) = ze3u_tf
225               e3v(ji,jj,jk,Nnn) = ze3v_tf
226               e3t(ji,jj,jk,Nnn) = ze3t_tf
227               !
228               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)   
229               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)     
230             END_3D   
231         ENDIF
232      ENDIF
233     
234      CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
235         &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
236
237!!an         
238
239      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
240      ! Set boundary conditions, time filter and swap time levels
241      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
242
243!!an TO BE ADDED : dyn_nxt
244!!                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
245!!an TO BE ADDED : a simplifier
246!                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
247 
248      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
249         !                                                  ! filtering "now" field
250         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) )
251      ENDIF
252     
253!!an
254
255
256      ! Swap time levels
257      Nrhs = Nbb
258      Nbb = Nnn
259      Nnn = Naa
260      Naa = Nrhs
261      !
262                         CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
263
264      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
265      ! diagnostics and outputs
266      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
267      IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats
268      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics
269   
270                        CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
271
272      !
273      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
274         
275
276#if defined key_agrif
277      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
278      ! AGRIF
279      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<     
280                         Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
281                         CALL Agrif_Integrate_ChildGrids( stp )       ! allows to finish all the Child Grids before updating
282
283                         IF( Agrif_NbStepint() == 0 ) THEN
284                            CALL Agrif_update_all( )                  ! Update all components
285                         ENDIF
286#endif
287      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
288
289      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
290      ! Control
291      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
292                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
293               
294                         
295      IF( kstp == nit000 ) THEN                          ! 1st time step only
296                                        CALL iom_close( numror )   ! close input  ocean restart file
297         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
298         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
299      ENDIF
300
301      !
302#if defined key_iomput
303      IF( kstp == nitend .OR. indic < 0 ) THEN
304                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
305                      IF(lrxios) CALL iom_context_finalize(      crxios_context          )
306      ENDIF
307#endif
308      !
309      IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
310         rDt = 2._wp * rn_Dt   
311         r1_Dt = 1._wp / rDt
312         l_1st_euler = .FALSE.     
313      ENDIF
314      !
315      IF( ln_timing )   CALL timing_stop('stp')
316      !
317   END SUBROUTINE stp
318   !
319   !!======================================================================
320END MODULE step
Note: See TracBrowser for help on using the repository browser.