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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3562

Last change on this file since 3562 was 3562, checked in by rblod, 11 years ago

Fix vvl and flux formulae, see ticket #952

  • Property svn:keywords set to Id
File size: 14.1 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   !!-------------------------------------------------------------------------
20 
21   !!-------------------------------------------------------------------------
22   !!   dyn_nxt      : obtain the next (after) horizontal velocity
23   !!-------------------------------------------------------------------------
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
26   USE sbc_oce         ! Surface boundary condition: ocean fields
27   USE phycst          ! physical constants
28   USE dynspg_oce      ! type of surface pressure gradient
29   USE dynadv          ! dynamics: vector invariant versus flux form
30   USE domvvl          ! variable volume
31   USE obc_oce         ! ocean open boundary conditions
32   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
33   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
34   USE obcvol          ! ocean open boundary condition (obc_vol routines)
35   USE bdy_oce         ! ocean open boundary conditions
36   USE bdydta          ! ocean open boundary conditions
37   USE bdydyn          ! ocean open boundary conditions
38   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
39   USE in_out_manager  ! I/O manager
40   USE lbclnk          ! lateral boundary condition (or mpp link)
41   USE lib_mpp         ! MPP library
42   USE wrk_nemo        ! Memory Allocation
43   USE prtctl          ! Print control
44#if defined key_agrif
45   USE agrif_opa_interp
46#endif
47   USE timing          ! Timing
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC    dyn_nxt   ! routine called by step.F90
53
54   !! * Substitutions
55#  include "domzgr_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE dyn_nxt ( kt )
64      !!----------------------------------------------------------------------
65      !!                  ***  ROUTINE dyn_nxt  ***
66      !!                   
67      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
68      !!             condition on the after velocity, achieved the time stepping
69      !!             by applying the Asselin filter on now fields and swapping
70      !!             the fields.
71      !!
72      !! ** Method  : * After velocity is compute using a leap-frog scheme:
73      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
74      !!             Note that with flux form advection and variable volume layer
75      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
76      !!             velocity.
77      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
78      !!             the time stepping has already been done in dynspg module
79      !!
80      !!              * Apply lateral boundary conditions on after velocity
81      !!             at the local domain boundaries through lbc_lnk call,
82      !!             at the one-way open boundaries (lk_obc=T),
83      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
84      !!
85      !!              * Apply the time filter applied and swap of the dynamics
86      !!             arrays to start the next time step:
87      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
88      !!                (un,vn) = (ua,va).
89      !!             Note that with flux form advection and variable volume layer
90      !!             (lk_vvl=T), the time filter is applied on thickness weighted
91      !!             velocity.
92      !!
93      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
94      !!               un,vn   now horizontal velocity of next time-step
95      !!----------------------------------------------------------------------
96      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
97      !
98      INTEGER  ::   ji, jj, jk   ! dummy loop indices
99      INTEGER  ::   iku, ikv     ! local integers
100#if ! defined key_dynspg_flt
101      REAL(wp) ::   z2dt         ! temporary scalar
102#endif
103      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
104      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
105      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
106      !!----------------------------------------------------------------------
107      !
108      IF( nn_timing == 1 )  CALL timing_start('dyn_nxt')
109      !
110      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
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      ! Next velocity :   Leap-frog time stepping
130      ! -------------
131      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
132      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
133      !
134      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
135         DO jk = 1, jpkm1
136            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
137            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
138         END DO
139      ELSE                                            ! applied on thickness weighted velocity
140         DO jk = 1, jpkm1
141            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
142               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
143               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
144            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
145               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
146               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
147         END DO
148      ENDIF
149
150
151      ! Update after velocity on domain lateral boundaries
152      ! --------------------------------------------------     
153      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
154      CALL lbc_lnk( va, 'V', -1. ) 
155      !
156# if defined key_obc
157      !                                !* OBC open boundaries
158      CALL obc_dyn( kt )
159      !
160      IF( .NOT. lk_dynspg_flt ) THEN
161         ! Flather boundary condition : - Update sea surface height on each open boundary
162         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
163         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
164         !                              - Correct the barotropic velocities
165         CALL obc_dyn_bt( kt )
166         !
167!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
168         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
169         !
170         IF( ln_vol_cst )   CALL obc_vol( kt )
171         !
172         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
173      ENDIF
174      !
175# elif defined key_bdy
176      !                                !* BDY open boundaries
177      IF( lk_dynspg_exp ) CALL bdy_dyn( kt )
178      IF( lk_dynspg_ts )  CALL bdy_dyn( kt, dyn3d_only=.true. )
179
180!!$   Do we need a call to bdy_vol here??
181      !
182# endif
183      !
184# if defined key_agrif
185      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
186# endif
187#endif
188
189      ! Time filter and swap of dynamics arrays
190      ! ------------------------------------------
191      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
192         DO jk = 1, jpkm1
193            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
194            vn(:,:,jk) = va(:,:,jk)
195         END DO
196      ELSE                                             !* Leap-Frog : Asselin filter and swap
197         !                                ! =============!
198         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
199            !                             ! =============!
200            DO jk = 1, jpkm1                             
201               DO jj = 1, jpj
202                  DO ji = 1, jpi   
203                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
204                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
205                     !
206                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
207                     vb(ji,jj,jk) = zvf
208                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
209                     vn(ji,jj,jk) = va(ji,jj,jk)
210                  END DO
211               END DO
212            END DO
213            !                             ! ================!
214         ELSE                             ! Variable volume !
215            !                             ! ================!
216            !
217            DO jk = 1, jpkm1                 ! Before scale factor at t-points
218               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
219                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
220                  &                         - 2._wp * fse3t_n(:,:,jk)            )
221            END DO
222            zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors
223            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
224            !
225            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation)
226               !
227               !                                      ! before scale factors at u- & v-pts (computed from fse3t_b)
228               CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) )
229               !
230               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity
231                  DO jj = 1, jpj                      !                                                 --------
232                     DO ji = 1, jpi
233                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
234                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
235                        !
236                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
237                        vb(ji,jj,jk) = zvf
238                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
239                        vn(ji,jj,jk) = va(ji,jj,jk)
240                     END DO
241                  END DO
242               END DO
243               !
244            ELSE                             ! flux form (thickness weighted calulation)
245               !
246               CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b)
247               !
248               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap:
249                  DO jj = 1, jpj                      !                   applied on thickness weighted velocity
250                     DO ji = 1, jpi                   !                              ---------------------------
251                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
252                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
253                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
254                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
255                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
256                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
257                        !
258                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
259                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
260                        !
261                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
262                        vb(ji,jj,jk) = zvf
263                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
264                        vn(ji,jj,jk) = va(ji,jj,jk)
265                     END DO
266                  END DO
267               END DO
268               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
269               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
270            ENDIF
271            !
272         ENDIF
273         !
274      ENDIF
275
276      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
277         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
278      !
279      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
280      !
281      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
282      !
283   END SUBROUTINE dyn_nxt
284
285   !!=========================================================================
286END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.