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

Last change on this file since 6004 was 6004, checked in by gm, 8 years ago

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