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

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

#1613: vvl by default: add ln_linssh and remove key_vvl
  • 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 non linear free surface,
76      !!             the leap-frog is applied on thickness weighted 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_bdy=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 non linear free surface,
90      !!             the time filter is applied on thickness weighted 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, z1_2dt   !   -      -
104      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
105      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
106      !!----------------------------------------------------------------------
107      !
108      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
109      !
110      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
111      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
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. ln_linssh ) 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) * e3u_b(:,:,jk)    &
145               &           + z2dt * ua(:,:,jk) * e3u_n(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
146            va(:,:,jk) = (          vb(:,:,jk) * e3v_b(:,:,jk)    &
147               &           + z2dt * va(:,:,jk) * e3v_n(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
148         END DO
149      ENDIF
150# endif
151
152# if defined key_dynspg_ts
153!!gm IF ( lk_dynspg_ts ) THEN ....
154      ! Ensure below that barotropic velocities match time splitting estimate
155      ! Compute actual transport and replace it with ts estimate at "after" time step
156      zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
157      zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
158      DO jk = 2, jpkm1
159         zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
160         zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
161      END DO
162      DO jk = 1, jpkm1
163         ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
164         va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
165      END DO
166
167      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
168         ! Remove advective velocity from "now velocities"
169         ! prior to asselin filtering     
170         ! In the forward case, this is done below after asselin filtering   
171         ! so that asselin contribution is removed at the same time
172         DO jk = 1, jpkm1
173            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
174            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
175         END DO 
176      ENDIF
177!!gm ENDIF
178# endif
179
180      ! Update after velocity on domain lateral boundaries
181      ! --------------------------------------------------     
182      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
183      CALL lbc_lnk( va, 'V', -1. ) 
184      !
185# if defined key_bdy
186      !                                !* BDY open boundaries
187      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
188      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
189
190!!$   Do we need a call to bdy_vol here??
191      !
192# endif
193      !
194# if defined key_agrif
195      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
196# endif
197#endif
198
199      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
200         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
201         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
202         !
203         !                                  ! Kinetic energy and Conversion
204         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
205         !
206         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
207            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
208            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
209            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
210            CALL iom_put( "vtrd_tot", zva )
211         ENDIF
212         !
213         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
214         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
215         !                                  !  computation of the asselin filter trends)
216      ENDIF
217
218      ! Time filter and swap of dynamics arrays
219      ! ------------------------------------------
220      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
221         DO jk = 1, jpkm1
222            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
223            vn(:,:,jk) = va(:,:,jk)
224         END DO
225         IF(.NOT.ln_linssh ) THEN
226            DO jk = 1, jpkm1
227               e3t_b(:,:,jk) = e3t_n(:,:,jk)
228               e3u_b(:,:,jk) = e3u_n(:,:,jk)
229               e3v_b(:,:,jk) = e3v_n(:,:,jk)
230            END DO
231         ENDIF
232      ELSE                                             !* Leap-Frog : Asselin filter and swap
233         !                                ! =============!
234         IF( ln_linssh ) THEN             ! Fixed volume !
235            !                             ! =============!
236            DO jk = 1, jpkm1                             
237               DO jj = 1, jpj
238                  DO ji = 1, jpi   
239                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
240                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
241                     !
242                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
243                     vb(ji,jj,jk) = zvf
244                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
245                     vn(ji,jj,jk) = va(ji,jj,jk)
246                  END DO
247               END DO
248            END DO
249            !                             ! ================!
250         ELSE                             ! Variable volume !
251            !                             ! ================!
252            ! Before scale factor at t-points
253            ! (used as a now filtered scale factor until the swap)
254            ! ----------------------------------------------------
255            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
256               ! No asselin filtering on thicknesses if forward time splitting
257                  e3t_b(:,:,:) = e3t_n(:,:,:)
258            ELSE
259               e3t_b(:,:,:) = e3t_n(:,:,:) + atfp * ( e3t_b(:,:,:) - 2._wp * e3t_n(:,:,:) + e3t_a(:,:,:) )
260               ! Add volume filter correction: compatibility with tracer advection scheme
261               ! => time filter + conservation correction (only at the first level)
262               IF ( nn_isf == 0) THEN   ! if no ice shelf melting
263                  e3t_b(:,:,1) = e3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
264                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
265               ELSE                     ! if ice shelf melting
266                  DO jj = 1,jpj
267                     DO ji = 1,jpi
268                        jk = mikt(ji,jj)
269                        e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       &
270                                          &                      * ( (emp_b(ji,jj)    - emp(ji,jj)   ) &
271                                          &                        - (rnf_b(ji,jj)    - rnf(ji,jj)   ) &
272                                          &                        + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk)
273                     END DO
274                  END DO
275               END IF
276            ENDIF
277            !
278            IF( ln_dynadv_vec ) THEN
279               ! Before scale factor at (u/v)-points
280               ! -----------------------------------
281               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
282               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
283               ! Leap-Frog - Asselin filter and swap: applied on velocity
284               ! -----------------------------------
285               DO jk = 1, jpkm1
286                  DO jj = 1, jpj
287                     DO ji = 1, jpi
288                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
289                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
290                        !
291                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
292                        vb(ji,jj,jk) = zvf
293                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
294                        vn(ji,jj,jk) = va(ji,jj,jk)
295                     END DO
296                  END DO
297               END DO
298               !
299            ELSE
300               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
301               !------------------------------------------------
302               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
303               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
304               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
305               ! -----------------------------------             ===========================
306               DO jk = 1, jpkm1
307                  DO jj = 1, jpj
308                     DO ji = 1, jpi                 
309                        zue3a = ua(ji,jj,jk) * e3u_a(ji,jj,jk)
310                        zve3a = va(ji,jj,jk) * e3v_a(ji,jj,jk)
311                        zue3n = un(ji,jj,jk) * e3u_n(ji,jj,jk)
312                        zve3n = vn(ji,jj,jk) * e3v_n(ji,jj,jk)
313                        zue3b = ub(ji,jj,jk) * e3u_b(ji,jj,jk)
314                        zve3b = vb(ji,jj,jk) * e3v_b(ji,jj,jk)
315                        !
316                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
317                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
318                        !
319                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
320                        vb(ji,jj,jk) = zvf
321                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
322                        vn(ji,jj,jk) = va(ji,jj,jk)
323                     END DO
324                  END DO
325               END DO
326               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
327               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
328            ENDIF
329            !
330         ENDIF
331         !
332         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
333            ! Revert "before" velocities to time split estimate
334            ! Doing it here also means that asselin filter contribution is removed 
335            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
336            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
337            DO jk = 2, jpkm1
338               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
339               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
340            END DO
341            DO jk = 1, jpkm1
342               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
343               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
344            END DO
345         ENDIF
346         !
347      ENDIF ! neuler =/0
348      !
349      ! Set "now" and "before" barotropic velocities for next time step:
350      ! JC: Would be more clever to swap variables than to make a full vertical
351      ! integration
352      !
353      !
354      IF(.NOT.ln_linssh ) THEN
355         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
356         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
357         DO jk = 2, jpkm1
358            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
359            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
360         END DO
361!!gm don't understand the use of umask_i ....
362         r1_hu_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
363         r1_hv_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
364      ENDIF
365      !
366      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
367      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
368      DO jk = 1, jpkm1
369         DO jj = 1, jpj
370            DO ji = 1, jpi
371               un_b(ji,jj) = un_b(ji,jj) + e3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
372               vn_b(ji,jj) = vn_b(ji,jj) + e3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
373               !
374               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
375               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
376            END DO
377         END DO
378      END DO
379      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
380      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
381      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
382      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
383      !
384      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
385         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
386         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
387         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
388      ENDIF
389      !
390      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
391         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
392      !
393      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
394      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
395      !
396      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
397      !
398   END SUBROUTINE dyn_nxt
399
400   !!=========================================================================
401END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.