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.
dynnxt.F90 in branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 4312

Last change on this file since 4312 was 4312, checked in by cetlod, 10 years ago

time-splitting corrections
  • Property svn:keywords set to Id
File size: 18.0 KB
Line 
1MODULE dynnxt
2   !!=========================================================================
3   !!                       ***  MODULE  dynnxt  ***
4   !! Ocean dynamics: time stepping
5   !!=========================================================================
6   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code
7   !!                 !  1990-10  (C. Levy, G. Madec)
8   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
9   !!            8.0  !  1997-02  (G. Madec & M. Imbard)  opa, release 8.0
10   !!            8.2  !  1997-04  (A. Weaver)  Euler forward step
11   !!             -   !  1997-06  (G. Madec)  lateral boudary cond., lbc routine
12   !!    NEMO    1.0  !  2002-08  (G. Madec)  F90: Free form and module
13   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond.
14   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
15   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.
16   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
17   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
18   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
19   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
20   !!-------------------------------------------------------------------------
21 
22   !!-------------------------------------------------------------------------
23   !!   dyn_nxt      : obtain the next (after) horizontal velocity
24   !!-------------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE sbc_oce         ! Surface boundary condition: ocean fields
28   USE phycst          ! physical constants
29   USE dynspg_oce      ! type of surface pressure gradient
30   USE dynadv          ! dynamics: vector invariant versus flux form
31   USE domvvl          ! variable volume
32   USE obc_oce         ! ocean open boundary conditions
33   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
34   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
35   USE obcvol          ! ocean open boundary condition (obc_vol routines)
36   USE bdy_oce         ! ocean open boundary conditions
37   USE bdydta          ! ocean open boundary conditions
38   USE bdydyn          ! ocean open boundary conditions
39   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
40   USE in_out_manager  ! I/O manager
41   USE lbclnk          ! lateral boundary condition (or mpp link)
42   USE lib_mpp         ! MPP library
43   USE wrk_nemo        ! Memory Allocation
44   USE prtctl          ! Print control
45   USE dynspg_ts       ! Barotropic velocities
46
47#if defined key_agrif
48   USE agrif_opa_interp
49#endif
50   USE timing          ! Timing
51
52   IMPLICIT NONE
53   PRIVATE
54
55   PUBLIC    dyn_nxt   ! routine called by step.F90
56
57   !! * Substitutions
58#  include "domzgr_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE dyn_nxt ( kt )
67      !!----------------------------------------------------------------------
68      !!                  ***  ROUTINE dyn_nxt  ***
69      !!                   
70      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
71      !!             condition on the after velocity, achieved the time stepping
72      !!             by applying the Asselin filter on now fields and swapping
73      !!             the fields.
74      !!
75      !! ** Method  : * After velocity is compute using a leap-frog scheme:
76      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
77      !!             Note that with flux form advection and variable volume layer
78      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
79      !!             velocity.
80      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
81      !!             the time stepping has already been done in dynspg module
82      !!
83      !!              * Apply lateral boundary conditions on after velocity
84      !!             at the local domain boundaries through lbc_lnk call,
85      !!             at the one-way open boundaries (lk_obc=T),
86      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
87      !!
88      !!              * Apply the time filter applied and swap of the dynamics
89      !!             arrays to start the next time step:
90      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
91      !!                (un,vn) = (ua,va).
92      !!             Note that with flux form advection and variable volume layer
93      !!             (lk_vvl=T), the time filter is applied on thickness weighted
94      !!             velocity.
95      !!
96      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
97      !!               un,vn   now horizontal velocity of next time-step
98      !!----------------------------------------------------------------------
99      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
100      !
101      INTEGER  ::   ji, jj, jk   ! dummy loop indices
102      INTEGER  ::   iku, ikv     ! local integers
103#if ! defined key_dynspg_flt
104      REAL(wp) ::   z2dt         ! temporary scalar
105#endif
106      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
107      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
108      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva, zhura, zhvra
109      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
110      !!----------------------------------------------------------------------
111      !
112      IF( nn_timing == 1 )  CALL timing_start('dyn_nxt')
113      !
114      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
115      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva, zhura, zhvra )
116      !
117      IF( kt == nit000 ) THEN
118         IF(lwp) WRITE(numout,*)
119         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
120         IF(lwp) WRITE(numout,*) '~~~~~~~'
121      ENDIF
122
123#if defined key_dynspg_flt
124      !
125      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
126      ! -------------
127
128      ! Update after velocity on domain lateral boundaries      (only local domain required)
129      ! --------------------------------------------------
130      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
131      CALL lbc_lnk( va, 'V', -1. ) 
132      !
133#else
134
135# if defined key_dynspg_exp
136      ! Next velocity :   Leap-frog time stepping
137      ! -------------
138      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
139      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
140      !
141      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
142         DO jk = 1, jpkm1
143            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
144            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
145         END DO
146      ELSE                                            ! applied on thickness weighted velocity
147         DO jk = 1, jpkm1
148            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
149               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
150               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
151            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
152               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
153               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
154         END DO
155      ENDIF
156# endif
157
158# if defined key_dynspg_ts
159      ! Ensure below that barotropic velocities match time splitting estimate
160      ! Compute actual transport and replace it with ts estimate at "after" time step
161      zua(:,:) = 0._wp
162      zva(:,:) = 0._wp
163      IF (lk_vvl) THEN
164         zhura(:,:) = 0._wp
165         zhvra(:,:) = 0._wp
166         DO jk = 1, jpkm1
167            zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
168            zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
169            zhura(:,:) = zhura(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
170            zhvra(:,:) = zhvra(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) 
171         END DO
172         zhura(:,:) = umask(:,:,1) / ( zhura(:,:) + 1._wp - umask(:,:,1) ) 
173         zhvra(:,:) = vmask(:,:,1) / ( zhvra(:,:) + 1._wp - vmask(:,:,1) ) 
174         DO jk = 1, jpkm1
175            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * zhura(:,:) + ua_b(:,:) ) * umask(:,:,jk)
176            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * zhvra(:,:) + va_b(:,:) ) * vmask(:,:,jk)
177         END DO
178      ELSE
179         DO jk = 1, jpkm1
180            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
181            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)   
182         END DO
183         DO jk = 1, jpkm1
184            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur(:,:) + ua_b(:,:) ) *umask(:,:,jk)
185            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr(:,:) + va_b(:,:) ) *vmask(:,:,jk)
186         END DO
187      ENDIF
188
189      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
190         ! Remove advective velocity from "now velocities"
191         ! prior to asselin filtering     
192         ! In the forward case, this is done below after asselin filtering   
193         ! so that asselin contribution is removed at the same time
194         DO jk = 1, jpkm1
195            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
196            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
197         END DO 
198      ENDIF
199# endif
200
201      ! Update after velocity on domain lateral boundaries
202      ! --------------------------------------------------     
203      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
204      CALL lbc_lnk( va, 'V', -1. ) 
205      !
206# if defined key_obc
207      !                                !* OBC open boundaries
208      IF( lk_obc ) CALL obc_dyn( kt )
209      !
210      IF( .NOT. lk_dynspg_flt ) THEN
211         ! Flather boundary condition : - Update sea surface height on each open boundary
212         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
213         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
214         !                              - Correct the barotropic velocities
215         IF( lk_obc ) CALL obc_dyn_bt( kt )
216         !
217!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
218         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
219         !
220         IF( lk_obc .AND. ln_vol_cst )   CALL obc_vol( kt )
221         !
222         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
223      ENDIF
224      !
225# elif defined key_bdy
226      !                                !* BDY open boundaries
227      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
228      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
229
230!!$   Do we need a call to bdy_vol here??
231      !
232# endif
233      !
234# if defined key_agrif
235      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
236# endif
237#endif
238
239      ! Time filter and swap of dynamics arrays
240      ! ------------------------------------------
241      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
242         DO jk = 1, jpkm1
243            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
244            vn(:,:,jk) = va(:,:,jk)
245         END DO
246         IF (lk_vvl) THEN
247            DO jk = 1, jpkm1
248               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
249               fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
250               fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
251            ENDDO
252         ENDIF
253      ELSE                                             !* Leap-Frog : Asselin filter and swap
254         !                                ! =============!
255         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
256            !                             ! =============!
257            DO jk = 1, jpkm1                             
258               DO jj = 1, jpj
259                  DO ji = 1, jpi   
260                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )
261                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )
262                     !
263                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
264                     vb(ji,jj,jk) = zvf
265                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
266                     vn(ji,jj,jk) = va(ji,jj,jk)
267                  END DO
268               END DO
269            END DO
270            !                             ! ================!
271         ELSE                             ! Variable volume !
272            !                             ! ================!
273            ! Before scale factor at t-points
274            ! (used as a now filtered scale factor until the swap)
275            ! ----------------------------------------------------
276            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
277               ! Remove asselin filtering on thicknesses if forward time splitting
278                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
279            ELSE
280               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
281               ! Add volume filter correction: compatibility with tracer advection scheme
282               ! => time filter + conservation correction (only at the first level)
283               fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
284            !
285            ENDIF
286            !
287            IF( ln_dynadv_vec ) THEN
288               ! Before scale factor at (u/v)-points
289               ! -----------------------------------
290               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
291               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
292               ! Leap-Frog - Asselin filter and swap: applied on velocity
293               ! -----------------------------------
294               DO jk = 1, jpkm1
295                  DO jj = 1, jpj
296                     DO ji = 1, jpi
297                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
298                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
299                        !
300                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
301                        vb(ji,jj,jk) = zvf
302                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
303                        vn(ji,jj,jk) = va(ji,jj,jk)
304                     END DO
305                  END DO
306               END DO
307               !
308            ELSE
309               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
310               !------------------------------------------------
311               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
312               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
313               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
314               ! -----------------------------------             ===========================
315               DO jk = 1, jpkm1
316                  DO jj = 1, jpj
317                     DO ji = 1, jpi                 
318                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
319                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
320                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
321                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
322                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
323                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
324                        !
325                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
326                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
327                        !
328                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
329                        vb(ji,jj,jk) = zvf
330                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
331                        vn(ji,jj,jk) = va(ji,jj,jk)
332                     END DO
333                  END DO
334               END DO
335               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
336               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
337            ENDIF
338            !
339         ENDIF
340         !
341         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
342            ! Revert "before" velocities to time split estimate
343            ! Doing it here also means that asselin filter contribution is removed 
344            zua(:,:) = 0._wp
345            zva(:,:) = 0._wp
346            IF (lk_vvl) THEN
347               DO jk = 1, jpkm1
348                  zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
349                  zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
350               END DO
351            ELSE
352               DO jk = 1, jpkm1
353                  zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
354                  zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
355               END DO
356            ENDIF
357            DO jk = 1, jpkm1
358               ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
359               vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
360            END DO
361         ENDIF
362         !
363      ENDIF ! neuler =/0
364
365      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
366         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
367      !
368      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
369      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva, zhura, zhvra )
370      !
371      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
372      !
373   END SUBROUTINE dyn_nxt
374
375   !!=========================================================================
376END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.