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_r12072_MERGE_OPTION2_2019/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DYN/dynnxt.F90 @ 12210

Last change on this file since 12210 was 12210, checked in by cetlod, 4 years ago

dev_merge_option2 : merge in fix_sn_cfctl_ticket2328 branch
  • Property svn:keywords set to Id
File size: 18.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.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 phycst         ! physical constants
32   USE dynadv         ! dynamics: vector invariant versus flux form
33   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
34   USE domvvl         ! variable volume
35   USE bdy_oce   , ONLY: ln_bdy
36   USE isf_oce   , ONLY: ln_isf     ! ice shelf
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   USE isfdynnxt , ONLY: isf_dynnxt ! ice shelf volume filter correction subroutine
44   !
45   USE in_out_manager ! I/O manager
46   USE iom            ! I/O manager library
47   USE lbclnk         ! lateral boundary condition (or mpp link)
48   USE lib_mpp        ! MPP library
49   USE prtctl         ! Print control
50   USE timing         ! Timing
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(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
103      !!----------------------------------------------------------------------
104      !
105      IF( ln_timing    )   CALL timing_start('dyn_nxt')
106      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
107      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
108      !
109      IF( kt == nit000 ) THEN
110         IF(lwp) WRITE(numout,*)
111         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
112         IF(lwp) WRITE(numout,*) '~~~~~~~'
113      ENDIF
114
115      IF ( ln_dynspg_ts ) THEN
116         ! Ensure below that barotropic velocities match time splitting estimate
117         ! Compute actual transport and replace it with ts estimate at "after" time step
118         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
119         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
120         DO jk = 2, jpkm1
121            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
122            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
123         END DO
124         DO jk = 1, jpkm1
125            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
126            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
127         END DO
128         !
129         IF( .NOT.ln_bt_fw ) THEN
130            ! Remove advective velocity from "now velocities"
131            ! prior to asselin filtering     
132            ! In the forward case, this is done below after asselin filtering   
133            ! so that asselin contribution is removed at the same time
134            DO jk = 1, jpkm1
135               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
136               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
137            END DO 
138         ENDIF
139      ENDIF
140
141      ! Update after velocity on domain lateral boundaries
142      ! --------------------------------------------------     
143# if defined key_agrif
144      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
145# endif
146      !
147      CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
148      !
149      !                                !* BDY open boundaries
150      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
151      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
152
153!!$   Do we need a call to bdy_vol here??
154      !
155      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
156         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
157         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
158         !
159         !                                  ! Kinetic energy and Conversion
160         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
161         !
162         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
163            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
164            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
165            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
166            CALL iom_put( "vtrd_tot", zva )
167         ENDIF
168         !
169         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
170         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
171         !                                  !  computation of the asselin filter trends)
172      ENDIF
173
174      ! Time filter and swap of dynamics arrays
175      ! ------------------------------------------
176      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
177         DO jk = 1, jpkm1
178            ub(:,:,jk) = un(:,:,jk)                         ! ub <-- un
179            vb(:,:,jk) = vn(:,:,jk)
180            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
181            vn(:,:,jk) = va(:,:,jk)
182         END DO
183         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
184!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a 
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               !
190               e3t_n(:,:,jk) = e3t_a(:,:,jk)
191               e3u_n(:,:,jk) = e3u_a(:,:,jk)
192               e3v_n(:,:,jk) = e3v_a(:,:,jk)
193            END DO
194!!gm BUG end
195         ENDIF
196                                                            !
197         
198      ELSE                                             !* Leap-Frog : Asselin filter and swap
199         !                                ! =============!
200         IF( ln_linssh ) THEN             ! Fixed volume !
201            !                             ! =============!
202            DO jk = 1, jpkm1                             
203               DO jj = 1, jpj
204                  DO ji = 1, jpi   
205                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
206                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
207                     !
208                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
209                     vb(ji,jj,jk) = zvf
210                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
211                     vn(ji,jj,jk) = va(ji,jj,jk)
212                  END DO
213               END DO
214            END DO
215            !                             ! ================!
216         ELSE                             ! Variable volume !
217            !                             ! ================!
218            ! Before scale factor at t-points
219            ! (used as a now filtered scale factor until the swap)
220            ! ----------------------------------------------------
221            DO jk = 1, jpkm1
222               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
223            END DO
224            ! Add volume filter correction: compatibility with tracer advection scheme
225            ! => time filter + conservation correction (only at the first level)
226            zcoef = atfp * rdt * r1_rau0
227
228            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
229
230            IF ( ln_rnf ) THEN
231               IF( ln_rnf_depth ) THEN
232                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
233                     DO jj = 1, jpj
234                        DO ji = 1, jpi
235                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
236                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
237                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
238                           ENDIF
239                        ENDDO
240                     ENDDO
241                  ENDDO
242               ELSE
243                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
244               ENDIF
245            END IF
246            !
247            ! ice shelf melting (deal separatly as it can be in depth)
248            ! PM: we could probably define a generic subroutine to do the in depth correction
249            !     to manage rnf, isf and possibly in the futur icb, tide water glacier (...)
250            !     ...(kt, coef, ktop, kbot, hz, fwf_b, fwf)
251            IF ( ln_isf ) CALL isf_dynnxt( kt, atfp * rdt )
252            !
253            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
254               ! Before filtered scale factor at (u/v)-points
255               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
256               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
257               DO jk = 1, jpkm1
258                  DO jj = 1, jpj
259                     DO ji = 1, jpi
260                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
261                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
262                        !
263                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
264                        vb(ji,jj,jk) = zvf
265                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
266                        vn(ji,jj,jk) = va(ji,jj,jk)
267                     END DO
268                  END DO
269               END DO
270               !
271            ELSE                          ! Asselin filter applied on thickness weighted velocity
272               !
273               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
274               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
275               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
276               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
277               DO jk = 1, jpkm1
278                  DO jj = 1, jpj
279                     DO ji = 1, jpi                 
280                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
281                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
282                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
283                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
284                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
285                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
286                        !
287                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
288                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
289                        !
290                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
291                        vb(ji,jj,jk) = zvf
292                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
293                        vn(ji,jj,jk) = va(ji,jj,jk)
294                     END DO
295                  END DO
296               END DO
297               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
298               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
299               !
300               DEALLOCATE( ze3u_f , ze3v_f )
301            ENDIF
302            !
303         ENDIF
304         !
305         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
306            ! Revert "before" velocities to time split estimate
307            ! Doing it here also means that asselin filter contribution is removed 
308            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
309            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
310            DO jk = 2, jpkm1
311               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
312               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
313            END DO
314            DO jk = 1, jpkm1
315               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
316               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
317            END DO
318         ENDIF
319         !
320      ENDIF ! neuler =/0
321      !
322      ! Set "now" and "before" barotropic velocities for next time step:
323      ! JC: Would be more clever to swap variables than to make a full vertical
324      ! integration
325      !
326      !
327      IF(.NOT.ln_linssh ) THEN
328         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
329         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
330         DO jk = 2, jpkm1
331            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
332            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
333         END DO
334         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
335         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
336      ENDIF
337      !
338      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
339      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
340      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
341      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
342      DO jk = 2, jpkm1
343         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
344         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
345         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
346         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
347      END DO
348      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
349      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
350      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
351      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
352      !
353      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
354         CALL iom_put(  "ubar", un_b(:,:) )
355         CALL iom_put(  "vbar", vn_b(:,:) )
356      ENDIF
357      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
358         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
359         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
360         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
361      ENDIF
362      !
363      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
364         &                                  tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
365      !
366      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
367      IF( l_trddyn     )   DEALLOCATE( zua, zva )
368      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
369      !
370   END SUBROUTINE dyn_nxt
371
372   !!=========================================================================
373END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.