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 NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/OCE/DYN/dynnxt.F90 @ 13662

Last change on this file since 13662 was 13662, checked in by clem, 3 years ago

update to almost r4.0.4
  • Property svn:keywords set to Id
File size: 19.3 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 sbcrnf         ! river runoffs
31   USE sbcisf         ! ice shelf
32   USE phycst         ! physical constants
33   USE dynadv         ! dynamics: vector invariant versus flux form
34   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
35   USE domvvl         ! variable volume
36   USE bdy_oce , ONLY : ln_bdy
37   USE bdydta         ! ocean open boundary conditions
38   USE bdydyn         ! ocean open boundary conditions
39   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
40   USE trd_oce        ! trends: ocean variables
41   USE trddyn         ! trend manager: dynamics
42   USE trdken         ! trend manager: kinetic energy
43   !
44   USE in_out_manager ! I/O manager
45   USE iom            ! I/O manager library
46   USE lbclnk         ! lateral boundary condition (or mpp link)
47   USE lib_mpp        ! MPP library
48   USE prtctl         ! Print control
49   USE timing         ! Timing
50   USE zdfdrg ,  ONLY : ln_drgice_imp, rCdU_top
51#if defined key_agrif
52   USE agrif_oce_interp
53#endif
54
55   IMPLICIT NONE
56   PRIVATE
57
58   PUBLIC    dyn_nxt   ! routine called by step.F90
59
60   !!----------------------------------------------------------------------
61   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
62   !! $Id$
63   !! Software governed by the CeCILL license (see ./LICENSE)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE dyn_nxt ( kt )
68      !!----------------------------------------------------------------------
69      !!                  ***  ROUTINE dyn_nxt  ***
70      !!                   
71      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
72      !!             condition on the after velocity, achieve the time stepping
73      !!             by applying the Asselin filter on now fields and swapping
74      !!             the fields.
75      !!
76      !! ** Method  : * Ensure after velocities transport matches time splitting
77      !!             estimate (ln_dynspg_ts=T)
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 (ln_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 non linear free surface,
89      !!             the time filter is applied on thickness weighted velocity.
90      !!             As a result, dyn_nxt MUST be called after tra_nxt.
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  ::   ikt          ! local integers
99      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
100      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
101      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
102      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zutau, zvtau
103      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
104      !!----------------------------------------------------------------------
105      !
106      IF( ln_timing    )   CALL timing_start('dyn_nxt')
107      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
108      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
109      !
110      IF( kt == nit000 ) THEN
111         IF(lwp) WRITE(numout,*)
112         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
113         IF(lwp) WRITE(numout,*) '~~~~~~~'
114      ENDIF
115
116      IF ( ln_dynspg_ts ) THEN
117         ! Ensure below that barotropic velocities match time splitting estimate
118         ! Compute actual transport and replace it with ts estimate at "after" time step
119         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
120         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
121         DO jk = 2, jpkm1
122            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
123            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
124         END DO
125         DO jk = 1, jpkm1
126            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
127            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
128         END DO
129         !
130         IF( .NOT.ln_bt_fw ) THEN
131            ! Remove advective velocity from "now velocities"
132            ! prior to asselin filtering     
133            ! In the forward case, this is done below after asselin filtering   
134            ! so that asselin contribution is removed at the same time
135            DO jk = 1, jpkm1
136               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
137               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
138            END DO 
139         ENDIF
140      ENDIF
141
142      ! Update after velocity on domain lateral boundaries
143      ! --------------------------------------------------     
144# if defined key_agrif
145      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
146# endif
147      !
148      CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
149      !
150      !                                !* BDY open boundaries
151      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
152      IF( ln_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      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
157         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
158         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
159         !
160         !                                  ! Kinetic energy and Conversion
161         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
162         !
163         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
164            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
165            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
166            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
167            CALL iom_put( "vtrd_tot", zva )
168         ENDIF
169         !
170         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
171         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
172         !                                  !  computation of the asselin filter trends)
173      ENDIF
174
175      ! Time filter and swap of dynamics arrays
176      ! ------------------------------------------
177      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
178         DO jk = 1, jpkm1
179            ub(:,:,jk) = un(:,:,jk)                         ! ub <-- un
180            vb(:,:,jk) = vn(:,:,jk)
181            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
182            vn(:,:,jk) = va(:,:,jk)
183         END DO
184         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
185!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a 
186            DO jk = 1, jpkm1
187!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
188!               e3u_b(:,:,jk) = e3u_n(:,:,jk)
189!               e3v_b(:,:,jk) = e3v_n(:,:,jk)
190               !
191               e3t_n(:,:,jk) = e3t_a(:,:,jk)
192               e3u_n(:,:,jk) = e3u_a(:,:,jk)
193               e3v_n(:,:,jk) = e3v_a(:,:,jk)
194            END DO
195!!gm BUG end
196         ENDIF
197                                                            !
198         
199      ELSE                                             !* Leap-Frog : Asselin filter and swap
200         !                                ! =============!
201         IF( ln_linssh ) THEN             ! Fixed volume !
202            !                             ! =============!
203            DO jk = 1, jpkm1                             
204               DO jj = 1, jpj
205                  DO ji = 1, jpi   
206                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
207                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
208                     !
209                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
210                     vb(ji,jj,jk) = zvf
211                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
212                     vn(ji,jj,jk) = va(ji,jj,jk)
213                  END DO
214               END DO
215            END DO
216            !                             ! ================!
217         ELSE                             ! Variable volume !
218            !                             ! ================!
219            ! Before scale factor at t-points
220            ! (used as a now filtered scale factor until the swap)
221            ! ----------------------------------------------------
222            DO jk = 1, jpkm1
223               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
224            END DO
225            ! Add volume filter correction: compatibility with tracer advection scheme
226            ! => time filter + conservation correction (only at the first level)
227            zcoef = atfp * rdt * r1_rau0
228
229            DO jk = 1, jpkm1
230               e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,jk) & 
231                             &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
232            END DO
233
234            IF ( ln_rnf ) THEN
235               DO jk = 1, jpkm1
236                  e3t_b(:,:,jk) = e3t_b(:,:,jk) + zcoef * ( rnf_b(:,:) - rnf(:,:) ) * tmask(:,:,jk) & 
237                                &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
238               END DO
239            ENDIF
240
241            IF ( ln_isf ) THEN
242               DO jk = 1, jpkm1
243                  e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( fwfisf_b(:,:) - fwfisf(:,:) ) * tmask(:,:,jk) & 
244                                &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
245               END DO
246            ENDIF
247            !
248            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
249               ! Before filtered scale factor at (u/v)-points
250               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
251               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
252               DO jk = 1, jpkm1
253                  DO jj = 1, jpj
254                     DO ji = 1, jpi
255                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
256                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
257                        !
258                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
259                        vb(ji,jj,jk) = zvf
260                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
261                        vn(ji,jj,jk) = va(ji,jj,jk)
262                     END DO
263                  END DO
264               END DO
265               !
266            ELSE                          ! Asselin filter applied on thickness weighted velocity
267               !
268               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
269               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
270               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
271               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
272               DO jk = 1, jpkm1
273                  DO jj = 1, jpj
274                     DO ji = 1, jpi                 
275                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
276                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
277                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
278                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
279                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
280                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
281                        !
282                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
283                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
284                        !
285                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
286                        vb(ji,jj,jk) = zvf
287                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
288                        vn(ji,jj,jk) = va(ji,jj,jk)
289                     END DO
290                  END DO
291               END DO
292               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
293               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
294               !
295               DEALLOCATE( ze3u_f , ze3v_f )
296            ENDIF
297            !
298         ENDIF
299         !
300         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
301            ! Revert "before" velocities to time split estimate
302            ! Doing it here also means that asselin filter contribution is removed 
303            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
304            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
305            DO jk = 2, jpkm1
306               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
307               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
308            END DO
309            DO jk = 1, jpkm1
310               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
311               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
312            END DO
313         ENDIF
314         !
315      ENDIF ! neuler =/0
316      !
317      ! Set "now" and "before" barotropic velocities for next time step:
318      ! JC: Would be more clever to swap variables than to make a full vertical
319      ! integration
320      !
321      !
322      IF(.NOT.ln_linssh ) THEN
323         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
324         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
325         DO jk = 2, jpkm1
326            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
327            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
328         END DO
329         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
330         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
331      ENDIF
332      !
333      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
334      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
335      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
336      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
337      DO jk = 2, jpkm1
338         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
339         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
340         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
341         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
342      END DO
343      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
344      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
345      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
346      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
347      !
348      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
349         CALL iom_put(  "ubar", un_b(:,:) )
350         CALL iom_put(  "vbar", vn_b(:,:) )
351      ENDIF
352      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
353         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
354         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
355         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
356      ENDIF
357      !
358      IF ( iom_use("utau") ) THEN
359         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
360            ALLOCATE(zutau(jpi,jpj)) 
361            DO jj = 2, jpjm1
362               DO ji = 2, jpim1
363                  jk = miku(ji,jj) 
364                  zutau(ji,jj) = utau(ji,jj) & 
365                  &  + 0.5_wp * rau0 * (rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) * ua(ji,jj,jk) 
366               END DO
367            END DO
368            CALL lbc_lnk( 'dynnxt' , zutau, 'U', -1.)
369            CALL iom_put(  "utau", zutau(:,:) )
370            DEALLOCATE(zutau)
371         ELSE
372            CALL iom_put(  "utau", utau(:,:) )
373         ENDIF
374      ENDIF
375      !
376      IF ( iom_use("vtau") ) THEN
377         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
378            ALLOCATE(zvtau(jpi,jpj))
379            DO jj = 2, jpjm1
380               DO ji = 2, jpim1
381                  jk = mikv(ji,jj)
382                  zvtau(ji,jj) = vtau(ji,jj) &
383                  &  + 0.5_wp * rau0 * (rCdU_top(ji,jj+1)+rCdU_top(ji,jj)) * va(ji,jj,jk)
384               END DO
385            END DO
386            CALL lbc_lnk( 'dynnxt' , zvtau, 'V', -1.)
387            CALL iom_put(  "vtau", zvtau(:,:) )
388            DEALLOCATE(zvtau)
389         ELSE
390            CALL iom_put(  "vtau", vtau(:,:) )
391         ENDIF
392      ENDIF
393      !
394      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
395         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
396      !
397      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
398      IF( l_trddyn     )   DEALLOCATE( zua, zva )
399      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
400      !
401   END SUBROUTINE dyn_nxt
402
403   !!=========================================================================
404END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.