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

source: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3970

Last change on this file since 3970 was 3970, checked in by cbricaud, 11 years ago

Time splitting update, see ticket #1079
  • Property svn:keywords set to Id
File size: 18.2 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
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 )
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         DO jk = 1, jpkm1
165            zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
166            zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)   
167         END DO
168         DO jk = 1, jpkm1
169            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) &
170            & / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) ) + ua_b(:,:) ) * umask(:,:,jk)
171            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) & 
172            & / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) ) + va_b(:,:) ) * vmask(:,:,jk)
173         END DO
174      ELSE
175         DO jk = 1, jpkm1
176            zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
177            zva(:,:) = zva(:,:) + fse3v(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)   
178         END DO
179         DO jk = 1, jpkm1
180            ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur(:,:) + ua_b(:,:) ) *umask(:,:,jk)
181            va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr(:,:) + va_b(:,:) ) *vmask(:,:,jk)
182         END DO
183      ENDIF
184
185      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
186         ! Remove advective velocity from "now velocities"
187         ! prior to asselin filtering     
188         ! In the forward case, this is done below after asselin filtering   
189         DO jk = 1, jpkm1
190            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
191            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
192         END DO 
193      ENDIF
194#endif
195
196      ! Update after velocity on domain lateral boundaries
197      ! --------------------------------------------------     
198      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
199      CALL lbc_lnk( va, 'V', -1. ) 
200      !
201# if defined key_obc
202      !                                !* OBC open boundaries
203      IF( lk_obc ) CALL obc_dyn( kt )
204      !
205      IF( .NOT. lk_dynspg_flt ) THEN
206         ! Flather boundary condition : - Update sea surface height on each open boundary
207         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
208         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
209         !                              - Correct the barotropic velocities
210         IF( lk_obc ) CALL obc_dyn_bt( kt )
211         !
212!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
213         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
214         !
215         IF( lk_obc .AND. ln_vol_cst )   CALL obc_vol( kt )
216         !
217         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
218      ENDIF
219      !
220# elif defined key_bdy
221      !                                !* BDY open boundaries
222      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
223      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
224
225!!$   Do we need a call to bdy_vol here??
226      !
227# endif
228      !
229# if defined key_agrif
230      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
231# endif
232#endif
233
234      ! Time filter and swap of dynamics arrays
235      ! ------------------------------------------
236      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
237         DO jk = 1, jpkm1
238            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
239            vn(:,:,jk) = va(:,:,jk)
240         END DO
241         IF (lk_vvl) THEN
242            DO jk = 1, jpkm1
243               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
244               fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
245               fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
246            ENDDO
247         ENDIF
248      ELSE                                             !* Leap-Frog : Asselin filter and swap
249         !                                ! =============!
250         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
251            !                             ! =============!
252            DO jk = 1, jpkm1                             
253               DO jj = 1, jpj
254                  DO ji = 1, jpi   
255                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
256                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
257                     !
258                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
259                     vb(ji,jj,jk) = zvf
260                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
261                     vn(ji,jj,jk) = va(ji,jj,jk)
262                  END DO
263               END DO
264            END DO
265            !                             ! ================!
266         ELSE                             ! Variable volume !
267            !                             ! ================!
268            !
269            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
270               ! Remove asselin filtering on thicknesses if forward time splitting
271               DO jk = 1, jpkm1   
272                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
273               ENDDO
274            ELSE
275
276               DO jk = 1, jpkm1                 ! Before scale factor at t-points
277                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
278                     &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
279                     &                         - 2._wp * fse3t_n(:,:,jk)            )
280               END DO
281               zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors
282               fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
283            ENDIF
284            !
285            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation)
286               !
287               !                                      ! before scale factors at u- & v-pts (computed from fse3t_b)
288               CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) )
289               !
290               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity
291                  DO jj = 1, jpj                      !                                                 --------
292                     DO ji = 1, jpi
293                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
294                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
295                        !
296                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
297                        vb(ji,jj,jk) = zvf
298                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
299                        vn(ji,jj,jk) = va(ji,jj,jk)
300                     END DO
301                  END DO
302               END DO
303               !
304            ELSE                             ! flux form (thickness weighted calulation)
305               !
306               CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b)
307               !
308               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap:
309                  DO jj = 1, jpj                      !                   applied on thickness weighted velocity
310                     DO ji = 1, jpi                   !                              ---------------------------
311                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
312                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
313                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
314                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
315                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
316                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
317                        !
318                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
319                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
320                        !
321                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
322                        vb(ji,jj,jk) = zvf
323                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
324                        vn(ji,jj,jk) = va(ji,jj,jk)
325                     END DO
326                  END DO
327               END DO
328               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
329               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
330            ENDIF
331            !
332         ENDIF
333         !
334#if defined key_dynspg_ts
335         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
336         ! Remove asselin filtering of barotropic velocities if forward time splitting
337         ! note that we replace barotropic velocities by advective velocities       
338            zua(:,:) = 0._wp
339            zva(:,:) = 0._wp
340            IF (lk_vvl) THEN
341               DO jk = 1, jpkm1
342                  zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
343                  zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
344               END DO
345               DO jk = 1, jpkm1
346                  ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) &
347                  & / ( hu_0(:,:) + sshu_n(:,:) + 1._wp - umask(:,:,1) ) - un_b(:,:)) * umask(:,:,jk)
348                  vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) & 
349                  & / ( hv_0(:,:) + sshv_n(:,:) + 1._wp - vmask(:,:,1) ) - vn_b(:,:)) * 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               DO jk = 1, jpkm1
357                  ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
358                  vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
359               END DO
360            ENDIF
361         ENDIF
362#endif
363         !
364      ENDIF ! neuler =/0
365
366      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
367         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
368      !
369      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
370      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
371      !
372      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
373      !
374   END SUBROUTINE dyn_nxt
375
376   !!=========================================================================
377END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.