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/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

#1260 : see the associated wiki page for explanation
  • Property svn:keywords set to Id
File size: 17.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 bdy_oce         ! ocean open boundary conditions
33   USE bdydta          ! ocean open boundary conditions
34   USE bdydyn          ! ocean open boundary conditions
35   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
36   USE in_out_manager  ! I/O manager
37   USE lbclnk          ! lateral boundary condition (or mpp link)
38   USE lib_mpp         ! MPP library
39   USE wrk_nemo        ! Memory Allocation
40   USE prtctl          ! Print control
41
42#if defined key_agrif
43   USE agrif_opa_interp
44#endif
45   USE timing          ! Timing
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC    dyn_nxt   ! routine called by step.F90
51
52   !! * Substitutions
53#  include "domzgr_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE dyn_nxt ( kt )
62      !!----------------------------------------------------------------------
63      !!                  ***  ROUTINE dyn_nxt  ***
64      !!                   
65      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
66      !!             condition on the after velocity, achieved the time stepping
67      !!             by applying the Asselin filter on now fields and swapping
68      !!             the fields.
69      !!
70      !! ** Method  : * After velocity is compute using a leap-frog scheme:
71      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
72      !!             Note that with flux form advection and variable volume layer
73      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
74      !!             velocity.
75      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
76      !!             the time stepping has already been done in dynspg module
77      !!
78      !!              * Apply lateral boundary conditions on after velocity
79      !!             at the local domain boundaries through lbc_lnk call,
80      !!             at the one-way open boundaries (lk_bdy=T),
81      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
82      !!
83      !!              * Apply the time filter applied and swap of the dynamics
84      !!             arrays to start the next time step:
85      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
86      !!                (un,vn) = (ua,va).
87      !!             Note that with flux form advection and variable volume layer
88      !!             (lk_vvl=T), the time filter is applied on thickness weighted
89      !!             velocity.
90      !!
91      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
92      !!               un,vn   now horizontal velocity of next time-step
93      !!----------------------------------------------------------------------
94      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
95      !
96      INTEGER  ::   ji, jj, jk   ! dummy loop indices
97      INTEGER  ::   iku, ikv     ! local integers
98#if ! defined key_dynspg_flt
99      REAL(wp) ::   z2dt         ! temporary scalar
100#endif
101      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
102      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
103      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva
104      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
105      !!----------------------------------------------------------------------
106      !
107      IF( nn_timing == 1 )  CALL timing_start('dyn_nxt')
108      !
109      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
110      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva )
111      !
112      IF( kt == nit000 ) THEN
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
115         IF(lwp) WRITE(numout,*) '~~~~~~~'
116      ENDIF
117
118#if defined key_dynspg_flt
119      !
120      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
121      ! -------------
122
123      ! Update after velocity on domain lateral boundaries      (only local domain required)
124      ! --------------------------------------------------
125      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
126      CALL lbc_lnk( va, 'V', -1. ) 
127      !
128#else
129
130# if defined key_dynspg_exp
131      ! Next velocity :   Leap-frog time stepping
132      ! -------------
133      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
134      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
135      !
136      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
137         DO jk = 1, jpkm1
138            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
139            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
140         END DO
141      ELSE                                            ! applied on thickness weighted velocity
142         DO jk = 1, jpkm1
143            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
144               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
145               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
146            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
147               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
148               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
149         END DO
150      ENDIF
151# endif
152
153# if defined key_dynspg_ts
154      ! Ensure below that barotropic velocities match time splitting estimate
155      ! Compute actual transport and replace it with ts estimate at "after" time step
156      zua(:,:) = 0._wp
157      zva(:,:) = 0._wp
158      DO jk = 1, jpkm1
159         zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
160         zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
161      END DO
162      DO jk = 1, jpkm1
163         ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
164         va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
165      END DO
166
167      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
168         ! Remove advective velocity from "now velocities"
169         ! prior to asselin filtering     
170         ! In the forward case, this is done below after asselin filtering   
171         ! so that asselin contribution is removed at the same time
172         DO jk = 1, jpkm1
173            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
174            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
175         END DO 
176      ENDIF
177# endif
178
179      ! Update after velocity on domain lateral boundaries
180      ! --------------------------------------------------     
181      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
182      CALL lbc_lnk( va, 'V', -1. ) 
183      !
184# if defined key_bdy
185      !                                !* BDY open boundaries
186      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
187      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
188
189!!$   Do we need a call to bdy_vol here??
190      !
191# endif
192      !
193# if defined key_agrif
194      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
195# endif
196#endif
197
198      ! Time filter and swap of dynamics arrays
199      ! ------------------------------------------
200      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
201         DO jk = 1, jpkm1
202            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
203            vn(:,:,jk) = va(:,:,jk)
204         END DO
205         IF (lk_vvl) THEN
206            DO jk = 1, jpkm1
207               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
208               fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
209               fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
210            ENDDO
211         ENDIF
212      ELSE                                             !* Leap-Frog : Asselin filter and swap
213         !                                ! =============!
214         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
215            !                             ! =============!
216            DO jk = 1, jpkm1                             
217               DO jj = 1, jpj
218                  DO ji = 1, jpi   
219                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )
220                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )
221                     !
222                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
223                     vb(ji,jj,jk) = zvf
224                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
225                     vn(ji,jj,jk) = va(ji,jj,jk)
226                  END DO
227               END DO
228            END DO
229            !                             ! ================!
230         ELSE                             ! Variable volume !
231            !                             ! ================!
232            ! Before scale factor at t-points
233            ! (used as a now filtered scale factor until the swap)
234            ! ----------------------------------------------------
235            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
236               ! No asselin filtering on thicknesses if forward time splitting
237                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
238            ELSE
239               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
240               ! Add volume filter correction: compatibility with tracer advection scheme
241               ! => time filter + conservation correction (only at the first level)
242               fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
243            !
244            ENDIF
245            !
246            IF( ln_dynadv_vec ) THEN
247               ! Before scale factor at (u/v)-points
248               ! -----------------------------------
249               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
250               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
251               ! Leap-Frog - Asselin filter and swap: applied on velocity
252               ! -----------------------------------
253               DO jk = 1, jpkm1
254                  DO jj = 1, jpj
255                     DO ji = 1, jpi
256                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
257                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
258                        !
259                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
260                        vb(ji,jj,jk) = zvf
261                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
262                        vn(ji,jj,jk) = va(ji,jj,jk)
263                     END DO
264                  END DO
265               END DO
266               !
267            ELSE
268               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
269               !------------------------------------------------
270               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
271               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
272               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
273               ! -----------------------------------             ===========================
274               DO jk = 1, jpkm1
275                  DO jj = 1, jpj
276                     DO ji = 1, jpi                 
277                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
278                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
279                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
280                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
281                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
282                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
283                        !
284                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
285                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
286                        !
287                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
288                        vb(ji,jj,jk) = zvf
289                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
290                        vn(ji,jj,jk) = va(ji,jj,jk)
291                     END DO
292                  END DO
293               END DO
294               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
295               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
296            ENDIF
297            !
298         ENDIF
299         !
300         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
301            ! Revert "before" velocities to time split estimate
302            ! Doing it here also means that asselin filter contribution is removed 
303            zua(:,:) = 0._wp
304            zva(:,:) = 0._wp
305            DO jk = 1, jpkm1
306               zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
307               zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
308            END DO
309            DO jk = 1, jpkm1
310               ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
311               vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
312            END DO
313         ENDIF
314         !
315      ENDIF ! neuler =/0
316      !
317      ! Set "now" and "before" barotropic velocities for next time step:
318      ! JC: Would be more clever to swap variables than to make a full vertical
319      ! integration
320      !
321      !
322      IF (lk_vvl) THEN
323         hu_b(:,:) = 0.
324         hv_b(:,:) = 0.
325         DO jk = 1, jpkm1
326            hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
327            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
328         END DO
329         hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )
330         hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )
331      ENDIF
332      !
333      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
334      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
335      !
336      DO jk = 1, jpkm1
337         DO jj = 1, jpj
338            DO ji = 1, jpi
339               un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
340               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
341               !
342               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
343               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
344            END DO
345         END DO
346      END DO
347      !
348      !
349      un_b(:,:) = un_b(:,:) * hur_a(:,:)
350      vn_b(:,:) = vn_b(:,:) * hvr_a(:,:)
351      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
352      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
353      !
354      !
355      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
356         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
357      !
358      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
359      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
360      !
361      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
362      !
363   END SUBROUTINE dyn_nxt
364
365   !!=========================================================================
366END MODULE dynnxt
367
Note: See TracBrowser for help on using the repository browser.