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

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

fixes from last merge qco r12983

File size: 17.5 KB
Line 
1MODULE stpLF
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   USE domqco
18
19   IMPLICIT NONE
20   PRIVATE
21
22   PUBLIC   stp_LF   ! 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_LF( )
41      INTEGER             ::   kstp   ! ocean time-step index
42#else
43   SUBROUTINE stp_LF( kstp )
44      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
45#endif
46      !!----------------------------------------------------------------------
47      !!                     ***  ROUTINE stp  ***
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     ! 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_LF')
82      !
83      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
84      ! model timestep
85      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
86      !
87      IF( l_1st_euler ) THEN     ! start or restart with Euler 1st time-step
88         rDt   =  rn_Dt   
89         r1_Dt = 1._wp / rDt
90      ENDIF
91
92      IF ( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity one for all to zero
93
94      !
95      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
96      ! update I/O and calendar
97      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
98                             indic = 0                ! reset to no error condition
99                             
100      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
101                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
102         IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
103                             CALL iom_init_closedef
104         IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
105      ENDIF
106      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
107                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
108      IF( ln_crs         )   CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell IOM we are at time step kstp
109
110      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
111      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
112      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
113      IF( ln_tide    )   CALL tide_update( kstp )                     ! update tide potential
114      IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
115      IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn )                   ! update dynamic & tracer data at open boundaries
116                         CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition
117
118      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
119      ! Ocean physics update
120      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
121      !  LATERAL  PHYSICS
122      !                                                                        ! eddy diffusivity coeff.
123      IF( l_ldfdyn_time )   CALL ldf_dyn( kstp, Nbb )                          ! eddy viscosity coeff.
124
125      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
126      !  Ocean dynamics : hdiv, ssh, e3, u, v, w
127      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
128
129                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
130                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
131                         vv(:,:,:,Nrhs) = 0._wp
132
133      IF( ln_bdy     )      CALL bdy_dyn3d_dmp ( kstp, Nbb,      uu, vv, Nrhs )  ! bdy damping trends
134
135#if defined key_agrif
136      IF(.NOT. Agrif_Root())  & 
137               &            CALL Agrif_Sponge_dyn        ! momentum sponge
138#endif
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 
143                            CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
144
145      IF( .NOT.ln_linssh )  CALL dom_qco_r3c   ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )   ! "after" ssh./h._0 ratio explicit
146      !IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt_st( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors
147!!an - calcul du gradient de pression horizontal (explicit)
148      DO_3D_00_00( 1, jpkm1 )
149         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
150         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
151      END_3D
152      !
153      ! add wind stress forcing and layer linear friction to the RHS
154      z1_2rho0 = 0.5_wp * r1_rho0
155      DO_3D_00_00(1,jpkm1)
156         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
157            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
158         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
159            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
160      END_3D   
161!!an         
162 
163      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
164      ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3
165      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
166!!st ssh_atf : add ssh filtering up there
167      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
168         !                                                  ! filtering "now" field
169         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2._wp * ssh(:,:,Nnn) + ssh(:,:,Naa) )
170      ENDIF
171!!st ssh_atf
172         
173!! what about  IF( .NOT.ln_linssh )  ?
174!!an futur module dyn_nxt (a la place de dyn_atf)
175     
176      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
177         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter)
178            DO_3D_00_00(1,jpkm1)
179               uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
180               vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
181             END_3D         
182         ELSE                             ! Leap Frog time stepping + Asselin filter         
183            DO_3D_11_11(1,jpkm1)
184               zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
185               zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
186               !                                                                  ! Asselin time filter on u,v (Nnn)
187               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)
188               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)
189               !             
190               uu(ji,jj,jk,Naa) = zua
191               vv(ji,jj,jk,Naa) = zva
192            END_3D
193            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn) )   ! "now" ssh/h_0 ratio from filtrered ssh
194#if ! defined key_qco
195            DO jk = 1, jpkm1
196               e3t(:,:,jk,Nnn) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Nnn) )
197               e3u(:,:,jk,Nnn) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Nnn) )
198               e3v(:,:,jk,Nnn) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Nnn) )
199            END DO
200#endif
201         ENDIF
202         !
203      ELSE                          ! flux form : applied on thickness weighted velocity
204         IF( l_1st_euler ) THEN           ! Euler time stepping (no Asselin filter)
205            DO_3D_00_00(1,jpkm1)
206               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
207               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
208               !                                                ! LF time stepping
209               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
210               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
211               !
212               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)   
213               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
214            END_3D
215         ELSE                             ! Leap Frog time stepping + Asselin filter
216            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:) )   ! "now" ssh/h_0 ratio from filtrered ssh
217            DO_3D_11_11(1,jpkm1)
218               zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn)
219               zve3n = e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nnn)
220               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
221               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
222               !                                                ! LF time stepping
223               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
224               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
225               !                                                ! Asselin time filter on e3u/v/t
226               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) ) 
227               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) )
228               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) ) 
229               !                                                ! Asselin time filter on u,v (Nnn)
230               uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / (e3u_0(ji,jj,jk) * ( 1._wp + r3u_f(ji,jj) ))
231               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / (e3v_0(ji,jj,jk) * ( 1._wp + r3v_f(ji,jj) ))
232               !
233               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)   
234               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
235            END_3D
236            r3t(:,:,Nnn) = r3t_f(:,:)
237            r3u(:,:,Nnn) = r3u_f(:,:)
238            r3v(:,:,Nnn) = r3v_f(:,:)
239#if ! defined key_qco
240            DO jk = 1, jpkm1
241               e3t(:,:,jk,Nnn) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Nnn) )
242               e3u(:,:,jk,Nnn) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Nnn) )
243               e3v(:,:,jk,Nnn) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Nnn) )
244            END DO
245#endif         
246         ENDIF
247      ENDIF
248
249
250      CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
251         &                       uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
252
253!!an         
254
255      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
256      ! Set boundary conditions, time filter and swap time levels
257      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
258
259!!an TO BE ADDED : dyn_nxt
260!!                         CALL dyn_atf       ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v  )  ! time filtering of "now" velocities and scale factors
261!!an TO BE ADDED : a simplifier
262!                         CALL ssh_atf       ( kstp, Nbb, Nnn, Naa, ssh )                     ! time filtering of "now" sea surface height
263!!st copied above
264!!      IF ( .NOT.( l_1st_euler ) ) THEN   ! Only do time filtering for leapfrog timesteps
265!!         !                                                  ! filtering "now" field
266!!         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) )
267!!      ENDIF
268!!st     
269!!an
270
271
272      ! Swap time levels
273      Nrhs = Nbb
274      Nbb = Nnn
275      Nnn = Naa
276      Naa = Nrhs
277      !
278!                         CALL dom_vvl_sf_update_st( kstp, Nbb, Nnn, Naa )  ! recompute vertical scale factors
279      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
280      ! diagnostics and outputs
281      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
282      IF( ln_floats  )   CALL flo_stp   ( kstp, Nbb, Nnn )      ! drifting Floats
283      IF( ln_diacfl  )   CALL dia_cfl   ( kstp,      Nnn )      ! Courant number diagnostics
284   
285                         CALL dia_wri   ( kstp,      Nnn )      ! ocean model: outputs
286
287      !
288      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
289         
290
291#if defined key_agrif
292      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
293      ! AGRIF
294      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<     
295                         Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
296                         CALL Agrif_Integrate_ChildGrids( stp_LF )       ! allows to finish all the Child Grids before updating
297
298                         IF( Agrif_NbStepint() == 0 ) THEN
299                            CALL Agrif_update_all( )                  ! Update all components
300                         ENDIF
301#endif
302      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
303
304      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
305      ! Control
306      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
307                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
308               
309                         
310      IF( kstp == nit000 ) THEN                          ! 1st time step only
311                                        CALL iom_close( numror )   ! close input  ocean restart file
312         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
313         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
314      ENDIF
315
316      !
317#if defined key_iomput
318      IF( kstp == nitend .OR. indic < 0 ) THEN
319                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
320                      IF(lrxios) CALL iom_context_finalize(      crxios_context          )
321      ENDIF
322#endif
323      !
324      IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
325         rDt = 2._wp * rn_Dt   
326         r1_Dt = 1._wp / rDt
327         l_1st_euler = .FALSE.     
328      ENDIF
329      !
330      IF( ln_timing )   CALL timing_stop('stp')
331      !
332   END SUBROUTINE stp_LF
333   !
334   !!======================================================================
335END MODULE stpLF
Note: See TracBrowser for help on using the repository browser.