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.
stpmlf.F90 in NEMO/trunk/src/SWE – NEMO

source: NEMO/trunk/src/SWE/stpmlf.F90 @ 14239

Last change on this file since 14239 was 14239, checked in by smasson, 3 years ago

trunk: replace key_iomput by key_xios
File size: 13.2 KB
Line 
1MODULE stpmlf
2   !!======================================================================
3   !!                       ***  MODULE stpMLF  ***
4   !! Time-stepping   : manager of the shallow water equation time stepping
5   !!                   Modified Leap Frog 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_MLF       : MLF Shallow Water Eq. time-stepping
13   !!----------------------------------------------------------------------
14   USE stp_oce        ! modules used in nemo_init and stp_MLF
15   !
16   USE domqco         ! quasi-eulerian coordinate
17   USE phycst         ! physical constants
18   USE usrdef_nam     ! user defined namelist parameters
19!!st   USE usrdef_sbc     ! user defined surface boundary cond
20   
21   IMPLICIT NONE
22   PRIVATE
23
24   PUBLIC   stp_MLF   ! called by nemogcm.F90
25   
26   !                                          !**  time level indices  **!
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   SUBROUTINE stp_MLF( kstp )
40      !!----------------------------------------------------------------------
41      !!                     ***  ROUTINE stp_MLF  ***
42      !!
43      !! ** Purpose : - Time stepping of shallow water (SHW) (momentum and ssh eqs.)
44      !!
45      !! ** Method  : -1- Update forcings
46      !!              -2- Update the ssh at Naa
47      !!              -3- Compute the momentum trends (Nrhs)
48      !!              -4- Update the horizontal velocity
49      !!              -5- Apply Asselin time filter to uu,vv,ssh
50      !!              -6- Outputs and diagnostics
51      !!----------------------------------------------------------------------
52      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
53      !
54      INTEGER ::   ji, jj, jk             ! dummy loop indices
55      INTEGER ::   indic                  ! error indicator if < 0
56      REAL(wp)::   z1_2rho0               ! local scalars
57      REAL(wp)::   zrhs_u, zue3a, zue3n, zue3b, zua   ! local scalars
58      REAL(wp)::   zrhs_v, zve3a, zve3n, zve3b, zva   !   -      -
59      !! ---------------------------------------------------------------------
60      !
61      IF( ln_timing )   CALL timing_start('stp_MLF')
62      !
63      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
64      ! model timestep
65      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
66      !
67      IF( l_1st_euler ) THEN   ! start or restart with Euler 1st time-step
68         rDt   = rn_Dt   
69         r1_Dt = 1._wp / rDt
70      ENDIF
71
72      IF( kstp == nit000 )   ww(:,:,:) = 0._wp   ! initialize vertical velocity once for all to zero
73
74      !
75      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
76      ! update I/O and calendar
77      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
78                             indic = 0                ! reset to no error condition
79
80      IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
81                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
82                             CALL iom_init_closedef
83      ENDIF
84      IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
85                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
86
87      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
88      ! Update external forcing   (SWE: surface boundary condition only)
89      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
90
91                             CALL sbc     ( kstp, Nbb, Nnn )                   ! Sea Boundary Condition
92
93      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
94      ! Ocean physics update   (SWE: eddy viscosity only)
95      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
96
97      IF( l_ldfdyn_time )   CALL ldf_dyn( kstp, Nbb )                          ! eddy viscosity coeff.
98
99      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
100      !  RHS of horizontal velocity
101      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
102
103                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
104                         vv(:,:,:,Nrhs) = 0._wp
105
106                         CALL dyn_adv( kstp, Nbb, Nnn, uu, vv, Nrhs )   ! advection (VF or FF)  ==> RHS
107                         CALL dyn_vor( kstp,      Nnn, uu, vv, Nrhs )   ! vorticity             ==> RHS
108                         CALL dyn_ldf( kstp, Nbb, Nnn, uu, vv, Nrhs )   ! lateral mixing                     ==> RHS
109
110      z1_2rho0 = 0.5_wp * r1_rho0
111      !
112      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
113         !                                                              ! horizontal pressure gradient
114         zrhs_u =        - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
115         zrhs_v =        - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
116         !                                                              ! wind stress and layer friction
117         zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn)   &
118            &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
119         zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn)   &
120            &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
121         !                                                              ! ==> RHS
122         uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
123         vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
124      END_3D
125
126      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
127      ! Time stepping of ssh Eq. (and update r3_Naa)
128      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
129      !        ! Leap Frog time stepping ==> ssh_Naa and r3_Naa
130                         CALL ssh_nxt( kstp, Nbb, Nnn, ssh   , Naa )    ! after ssh
131      !                                                                 ! after ssh/h_0 ratio explicit
132                         CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
133      !        ! Asselin filter ==> ssh_Nnn filtered
134      IF ( .NOT.( l_1st_euler ) ) THEN   ! Time filtering of now ssh
135         ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2._wp * ssh(:,:,Nnn) + ssh(:,:,Naa) )
136      ENDIF
137
138      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
139      ! Time stepping of dynamics (u,v)
140      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
141
142      IF( ln_dynadv_vec ) THEN      ! vector invariant form : applied on velocity
143         IF( l_1st_euler ) THEN          ! Euler time stepping (no Asselin filter)
144            DO_3D( 0,0, 0,0, 1,jpkm1)
145               uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
146               vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
147             END_3D         
148         ELSE                            ! Leap Frog time stepping + Asselin filter         
149            DO_3D( 0,0, 0,0, 1,jpkm1)
150               zua = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
151               zva = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
152               !                                                                  ! Asselin time filter on u,v (Nnn)
153               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)
154               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)
155               !             
156               uu(ji,jj,jk,Naa) = zua
157               vv(ji,jj,jk,Naa) = zva
158            END_3D
159            !                            ! Update r3_Nnn
160            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn) )   ! now ssh/h_0 ratio from filtrered ssh
161         ENDIF
162         !
163      ELSE                          ! flux form : applied on thickness weighted velocity
164         IF( l_1st_euler ) THEN          ! Euler time stepping (no Asselin filter)
165            DO_3D( 0,0, 0,0, 1,jpkm1)
166               zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
167               zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
168               !                                                ! Euler time stepping
169               zue3a = zue3b + rDt * e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
170               zve3a = zve3b + rDt * e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
171               !
172               uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)   
173               vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
174            END_3D
175         ELSE                             ! Leap Frog time stepping + Asselin filter
176            CALL dom_qco_r3c( ssh(:,:,Nnn), r3t_f(:,:), r3u_f(:,:), r3v_f(:,:) )   ! now ssh/h_0 ratio from filtrered ssh
177            DO_3D( 0,0, 0,0, 1,jpkm1)
178               zue3n = ( 1._wp + r3u(ji,jj,Nnn) ) * uu(ji,jj,jk,Nnn)
179               zve3n = ( 1._wp + r3v(ji,jj,Nnn) ) * vv(ji,jj,jk,Nnn)
180               zue3b = ( 1._wp + r3u(ji,jj,Nbb) ) * uu(ji,jj,jk,Nbb)
181               zve3b = ( 1._wp + r3v(ji,jj,Nbb) ) * vv(ji,jj,jk,Nbb)
182               !                                                ! LF time stepping
183               zue3a = zue3b + rDt * ( 1._wp + r3u(ji,jj,Nnn) ) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
184               zve3a = zve3b + rDt * ( 1._wp + r3v(ji,jj,Nnn) ) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
185               !                                                ! Asselin time filter on u,v (Nnn)
186               uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ( 1._wp + r3u_f(ji,jj) )
187               vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ( 1._wp + r3v_f(ji,jj) )
188               !
189               uu(ji,jj,jk,Naa) = zue3a / ( 1._wp + r3u(ji,jj,Naa) )   
190               vv(ji,jj,jk,Naa) = zve3a / ( 1._wp + r3v(ji,jj,Naa) )
191            END_3D
192            !                             ! Update r3_Nnn with time filtered values
193            r3t(:,:,Nnn) = r3t_f(:,:)
194            r3u(:,:,Nnn) = r3u_f(:,:)
195            r3v(:,:,Nnn) = r3v_f(:,:)
196         ENDIF
197      ENDIF
198
199      CALL lbc_lnk_multi( 'stp_MLF', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1.,   &   !* local domain boundaries
200         &                           uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1.    )     
201
202      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
203      ! Set boundary conditions, time filter and swap time levels
204      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
205
206      ! Swap time levels
207      Nrhs = Nbb
208      Nbb = Nnn
209      Nnn = Naa
210      Naa = Nrhs
211
212      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
213      ! diagnostics and outputs
214      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
215     
216      IF( ln_diacfl  )   CALL dia_cfl  ( kstp,      Nnn )      ! Courant number diagnostics
217                         CALL dia_wri  ( kstp,      Nnn )      ! ocean model: outputs
218      !
219      IF( lrst_oce   )   CALL rst_write( kstp, Nbb, Nnn )      ! write output ocean restart file
220
221      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
222      ! Control
223      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
224                         CALL stp_ctl_SWE   ( kstp, Nnn )
225
226      IF( kstp == nit000 ) THEN                          ! 1st time step only
227                                        CALL iom_close( numror )   ! close input  ocean restart file
228         IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
229         IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
230      ENDIF
231
232      !
233#if defined key_xios
234      IF( kstp == nitend .OR. indic < 0 ) THEN
235!!st : cxios_context needed ? because opened earlier ???         
236         CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
237      ENDIF
238#endif
239      !
240      IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
241         rDt   = 2._wp * rn_Dt   
242         r1_Dt = 1._wp / rDt
243         l_1st_euler = .FALSE.     
244      ENDIF
245      !
246      IF( ln_timing )   CALL timing_stop('stp_MLF')
247      !
248   END SUBROUTINE stp_MLF
249
250   !!======================================================================
251END MODULE stpmlf
Note: See TracBrowser for help on using the repository browser.