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 @ 4328

Last change on this file since 4328 was 4328, checked in by davestorkey, 10 years ago

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