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

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 5845

Last change on this file since 5845 was 5845, checked in by gm, 9 years ago

#1613: vvl by default: suppression of domzgr_substitute.h90

  • Property svn:keywords set to Id
File size: 19.5 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   !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
21   !!-------------------------------------------------------------------------
22 
23   !!-------------------------------------------------------------------------
24   !!   dyn_nxt      : obtain the next (after) horizontal velocity
25   !!-------------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE sbc_oce         ! Surface boundary condition: ocean fields
29   USE phycst          ! physical constants
30   USE dynspg_oce      ! type of surface pressure gradient
31   USE dynadv          ! dynamics: vector invariant versus flux form
32   USE domvvl          ! variable volume
33   USE bdy_oce         ! ocean open boundary conditions
34   USE bdydta          ! ocean open boundary conditions
35   USE bdydyn          ! ocean open boundary conditions
36   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
37   USE trd_oce         ! trends: ocean variables
38   USE trddyn          ! trend manager: dynamics
39   USE trdken          ! trend manager: kinetic energy
40   !
41   USE in_out_manager  ! I/O manager
42   USE iom             ! I/O manager library
43   USE lbclnk          ! lateral boundary condition (or mpp link)
44   USE lib_mpp         ! MPP library
45   USE wrk_nemo        ! Memory Allocation
46   USE prtctl          ! Print control
47   USE timing          ! Timing
48#if defined key_agrif
49   USE agrif_opa_interp
50#endif
51
52   IMPLICIT NONE
53   PRIVATE
54
55   PUBLIC    dyn_nxt   ! routine called by step.F90
56
57   !!----------------------------------------------------------------------
58   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
59   !! $Id$
60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64   SUBROUTINE dyn_nxt ( kt )
65      !!----------------------------------------------------------------------
66      !!                  ***  ROUTINE dyn_nxt  ***
67      !!                   
68      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
69      !!             condition on the after velocity, achieved the time stepping
70      !!             by applying the Asselin filter on now fields and swapping
71      !!             the fields.
72      !!
73      !! ** Method  : * After velocity is compute using a leap-frog scheme:
74      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
75      !!             Note that with flux form advection and variable volume layer
76      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
77      !!             velocity.
78      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
79      !!             the time stepping has already been done in dynspg module
80      !!
81      !!              * Apply lateral boundary conditions on after velocity
82      !!             at the local domain boundaries through lbc_lnk call,
83      !!             at the one-way open boundaries (lk_bdy=T),
84      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
85      !!
86      !!              * Apply the time filter applied and swap of the dynamics
87      !!             arrays to start the next time step:
88      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
89      !!                (un,vn) = (ua,va).
90      !!             Note that with flux form advection and variable volume layer
91      !!             (lk_vvl=T), the time filter is applied on thickness weighted
92      !!             velocity.
93      !!
94      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
95      !!               un,vn   now horizontal velocity of next time-step
96      !!----------------------------------------------------------------------
97      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
98      !
99      INTEGER  ::   ji, jj, jk   ! dummy loop indices
100      INTEGER  ::   iku, ikv     ! local integers
101#if ! defined key_dynspg_flt
102      REAL(wp) ::   z2dt         ! temporary scalar
103#endif
104      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars
105      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
106      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
107      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
108      !!----------------------------------------------------------------------
109      !
110      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
111      !
112      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
113      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
114      !
115      IF( kt == nit000 ) THEN
116         IF(lwp) WRITE(numout,*)
117         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
118         IF(lwp) WRITE(numout,*) '~~~~~~~'
119      ENDIF
120
121#if defined key_dynspg_flt
122      !
123      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
124      ! -------------
125
126      ! Update after velocity on domain lateral boundaries      (only local domain required)
127      ! --------------------------------------------------
128      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
129      CALL lbc_lnk( va, 'V', -1. ) 
130      !
131#else
132
133# if defined key_dynspg_exp
134      ! Next velocity :   Leap-frog time stepping
135      ! -------------
136      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
137      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
138      !
139      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
140         DO jk = 1, jpkm1
141            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
142            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
143         END DO
144      ELSE                                            ! applied on thickness weighted velocity
145         DO jk = 1, jpkm1
146            ua(:,:,jk) = (          ub(:,:,jk) * e3u_b(:,:,jk)    &
147               &           + z2dt * ua(:,:,jk) * e3u_n(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
148            va(:,:,jk) = (          vb(:,:,jk) * e3v_b(:,:,jk)    &
149               &           + z2dt * va(:,:,jk) * e3v_n(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
150         END DO
151      ENDIF
152# endif
153
154# if defined key_dynspg_ts
155!!gm IF ( lk_dynspg_ts ) THEN ....
156      ! Ensure below that barotropic velocities match time splitting estimate
157      ! Compute actual transport and replace it with ts estimate at "after" time step
158      zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
159      zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
160      DO jk = 2, jpkm1
161         zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
162         zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
163      END DO
164      DO jk = 1, jpkm1
165         ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
166         va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
167      END DO
168
169      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
170         ! Remove advective velocity from "now velocities"
171         ! prior to asselin filtering     
172         ! In the forward case, this is done below after asselin filtering   
173         ! so that asselin contribution is removed at the same time
174         DO jk = 1, jpkm1
175            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
176            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
177         END DO 
178      ENDIF
179!!gm ENDIF
180# endif
181
182      ! Update after velocity on domain lateral boundaries
183      ! --------------------------------------------------     
184      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
185      CALL lbc_lnk( va, 'V', -1. ) 
186      !
187# if defined key_bdy
188      !                                !* BDY open boundaries
189      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
190      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
191
192!!$   Do we need a call to bdy_vol here??
193      !
194# endif
195      !
196# if defined key_agrif
197      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
198# endif
199#endif
200
201      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
202         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
203         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
204         !
205         !                                  ! Kinetic energy and Conversion
206         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
207         !
208         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
209            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
210            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
211            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
212            CALL iom_put( "vtrd_tot", zva )
213         ENDIF
214         !
215         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
216         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
217         !                                  !  computation of the asselin filter trends)
218      ENDIF
219
220      ! Time filter and swap of dynamics arrays
221      ! ------------------------------------------
222      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
223         DO jk = 1, jpkm1
224            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
225            vn(:,:,jk) = va(:,:,jk)
226         END DO
227         IF (lk_vvl) THEN
228            DO jk = 1, jpkm1
229               e3t_b(:,:,jk) = e3t_n(:,:,jk)
230               e3u_b(:,:,jk) = e3u_n(:,:,jk)
231               e3v_b(:,:,jk) = e3v_n(:,:,jk)
232            ENDDO
233         ENDIF
234      ELSE                                             !* Leap-Frog : Asselin filter and swap
235         !                                ! =============!
236         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
237            !                             ! =============!
238            DO jk = 1, jpkm1                             
239               DO jj = 1, jpj
240                  DO ji = 1, jpi   
241                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
242                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
243                     !
244                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
245                     vb(ji,jj,jk) = zvf
246                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
247                     vn(ji,jj,jk) = va(ji,jj,jk)
248                  END DO
249               END DO
250            END DO
251            !                             ! ================!
252         ELSE                             ! Variable volume !
253            !                             ! ================!
254            ! Before scale factor at t-points
255            ! (used as a now filtered scale factor until the swap)
256            ! ----------------------------------------------------
257            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
258               ! No asselin filtering on thicknesses if forward time splitting
259                  e3t_b(:,:,:) = e3t_n(:,:,:)
260            ELSE
261               e3t_b(:,:,:) = e3t_n(:,:,:) + atfp * ( e3t_b(:,:,:) - 2._wp * e3t_n(:,:,:) + e3t_a(:,:,:) )
262               ! Add volume filter correction: compatibility with tracer advection scheme
263               ! => time filter + conservation correction (only at the first level)
264               IF ( nn_isf == 0) THEN   ! if no ice shelf melting
265                  e3t_b(:,:,1) = e3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
266                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
267               ELSE                     ! if ice shelf melting
268                  DO jj = 1,jpj
269                     DO ji = 1,jpi
270                        jk = mikt(ji,jj)
271                        e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       &
272                                          &                      * ( (emp_b(ji,jj)    - emp(ji,jj)   ) &
273                                          &                        - (rnf_b(ji,jj)    - rnf(ji,jj)   ) &
274                                          &                        + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk)
275                     END DO
276                  END DO
277               END IF
278            ENDIF
279            !
280            IF( ln_dynadv_vec ) THEN
281               ! Before scale factor at (u/v)-points
282               ! -----------------------------------
283               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
284               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
285               ! Leap-Frog - Asselin filter and swap: applied on velocity
286               ! -----------------------------------
287               DO jk = 1, jpkm1
288                  DO jj = 1, jpj
289                     DO ji = 1, jpi
290                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
291                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
292                        !
293                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
294                        vb(ji,jj,jk) = zvf
295                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
296                        vn(ji,jj,jk) = va(ji,jj,jk)
297                     END DO
298                  END DO
299               END DO
300               !
301            ELSE
302               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
303               !------------------------------------------------
304               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
305               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
306               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
307               ! -----------------------------------             ===========================
308               DO jk = 1, jpkm1
309                  DO jj = 1, jpj
310                     DO ji = 1, jpi                 
311                        zue3a = ua(ji,jj,jk) * e3u_a(ji,jj,jk)
312                        zve3a = va(ji,jj,jk) * e3v_a(ji,jj,jk)
313                        zue3n = un(ji,jj,jk) * e3u_n(ji,jj,jk)
314                        zve3n = vn(ji,jj,jk) * e3v_n(ji,jj,jk)
315                        zue3b = ub(ji,jj,jk) * e3u_b(ji,jj,jk)
316                        zve3b = vb(ji,jj,jk) * e3v_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               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
329               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
330            ENDIF
331            !
332         ENDIF
333         !
334         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
335            ! Revert "before" velocities to time split estimate
336            ! Doing it here also means that asselin filter contribution is removed 
337            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
338            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
339            DO jk = 2, jpkm1
340               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
341               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
342            END DO
343            DO jk = 1, jpkm1
344               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
345               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
346            END DO
347         ENDIF
348         !
349      ENDIF ! neuler =/0
350      !
351      ! Set "now" and "before" barotropic velocities for next time step:
352      ! JC: Would be more clever to swap variables than to make a full vertical
353      ! integration
354      !
355      !
356      IF (lk_vvl) THEN
357         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
358         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
359         DO jk = 2, jpkm1
360            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
361            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
362         END DO
363!!gm don't understand the use of umask_i ....
364         r1_hu_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
365         r1_hv_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
366      ENDIF
367      !
368      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
369      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
370      DO jk = 1, jpkm1
371         DO jj = 1, jpj
372            DO ji = 1, jpi
373               un_b(ji,jj) = un_b(ji,jj) + e3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
374               vn_b(ji,jj) = vn_b(ji,jj) + e3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
375               !
376               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
377               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
378            END DO
379         END DO
380      END DO
381      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
382      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
383      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
384      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
385      !
386      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
387         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
388         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
389         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
390      ENDIF
391      !
392      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
393         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
394      !
395      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
396      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
397      !
398      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
399      !
400   END SUBROUTINE dyn_nxt
401
402   !!=========================================================================
403END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.