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/2017/dev_merge_2017/NEMOGCM/NEMO/OCE_SRC/DYN – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OCE_SRC/DYN/dynnxt.F90 @ 9570

Last change on this file since 9570 was 9570, checked in by nicolasmartin, 6 years ago

Global renaming for core routines (./NEMO)

  • Folders
    • LIM_SRC_3 -> ICE_SRC
    • OPA_SRC -> OCE_SRC
  • CPP key: key_lim3 -> key_si3
  • Modules, (sub)routines and variables names
    • MPI: mpi_comm_opa -> mpi_comm_oce, MPI_COMM_OPA -> MPI_COMM_OCE, mpi_init_opa -> mpi_init_oce
    • AGRIF: agrif_opa_* -> agrif_oce_*, agrif_lim3_* -> agrif_si3_* and few more
    • TOP-PISCES: p.zlim -> p.zice, namp.zlim -> namp.zice
  • Comments
    • NEMO/OPA -> NEMO/OCE
    • ESIM|LIM3 -> SI3
  • 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 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#if defined key_agrif
51   USE agrif_oce_interp
52#endif
53
54   IMPLICIT NONE
55   PRIVATE
56
57   PUBLIC    dyn_nxt   ! routine called by step.F90
58
59   !!----------------------------------------------------------------------
60   !! NEMO/OCE 4.0 , NEMO Consortium (2017)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE dyn_nxt ( kt )
67      !!----------------------------------------------------------------------
68      !!                  ***  ROUTINE dyn_nxt  ***
69      !!                   
70      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
71      !!             condition on the after velocity, achieve the time stepping
72      !!             by applying the Asselin filter on now fields and swapping
73      !!             the fields.
74      !!
75      !! ** Method  : * Ensure after velocities transport matches time splitting
76      !!             estimate (ln_dynspg_ts=T)
77      !!
78      !!              * Apply lateral boundary conditions on after velocity
79      !!             at the local domain boundaries through lbc_lnk call,
80      !!             at the one-way open boundaries (ln_bdy=T),
81      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
82      !!
83      !!              * Apply the time filter applied and swap of the dynamics
84      !!             arrays to start the next time step:
85      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
86      !!                (un,vn) = (ua,va).
87      !!             Note that with flux form advection and non linear free surface,
88      !!             the time filter is applied on thickness weighted velocity.
89      !!             As a result, dyn_nxt MUST be called after tra_nxt.
90      !!
91      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
92      !!               un,vn   now horizontal velocity of next time-step
93      !!----------------------------------------------------------------------
94      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
95      !
96      INTEGER  ::   ji, jj, jk   ! dummy loop indices
97      INTEGER  ::   ikt          ! local integers
98      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
99      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
100      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
101      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
102      !!----------------------------------------------------------------------
103      !
104      IF( ln_timing    )   CALL timing_start('dyn_nxt')
105      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
106      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
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(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
135               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + 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_multi( ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
147      !
148      !                                !* BDY open boundaries
149      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
150      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
151
152!!$   Do we need a call to bdy_vol here??
153      !
154      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
155         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
156         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
157         !
158         !                                  ! Kinetic energy and Conversion
159         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
160         !
161         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
162            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
163            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
164            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
165            CALL iom_put( "vtrd_tot", zva )
166         ENDIF
167         !
168         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
169         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
170         !                                  !  computation of the asselin filter trends)
171      ENDIF
172
173      ! Time filter and swap of dynamics arrays
174      ! ------------------------------------------
175      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
176         DO jk = 1, jpkm1
177            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
178            vn(:,:,jk) = va(:,:,jk)
179         END DO
180         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
181!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a 
182            DO jk = 1, jpkm1
183!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
184!               e3u_b(:,:,jk) = e3u_n(:,:,jk)
185!               e3v_b(:,:,jk) = e3v_n(:,:,jk)
186               !
187               e3t_n(:,:,jk) = e3t_a(:,:,jk)
188               e3u_n(:,:,jk) = e3u_a(:,:,jk)
189               e3v_n(:,:,jk) = e3v_a(:,:,jk)
190            END DO
191!!gm BUG end
192         ENDIF
193                                                            !
194         
195      ELSE                                             !* Leap-Frog : Asselin filter and swap
196         !                                ! =============!
197         IF( ln_linssh ) THEN             ! Fixed volume !
198            !                             ! =============!
199            DO jk = 1, jpkm1                             
200               DO jj = 1, jpj
201                  DO ji = 1, jpi   
202                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
203                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
204                     !
205                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
206                     vb(ji,jj,jk) = zvf
207                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
208                     vn(ji,jj,jk) = va(ji,jj,jk)
209                  END DO
210               END DO
211            END DO
212            !                             ! ================!
213         ELSE                             ! Variable volume !
214            !                             ! ================!
215            ! Before scale factor at t-points
216            ! (used as a now filtered scale factor until the swap)
217            ! ----------------------------------------------------
218            DO jk = 1, jpkm1
219               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
220            END DO
221            ! Add volume filter correction: compatibility with tracer advection scheme
222            ! => time filter + conservation correction (only at the first level)
223            zcoef = atfp * rdt * r1_rau0
224
225            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
226
227            IF ( ln_rnf ) THEN
228               IF( ln_rnf_depth ) THEN
229                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
230                     DO jj = 1, jpj
231                        DO ji = 1, jpi
232                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
233                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
234                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
235                           ENDIF
236                        ENDDO
237                     ENDDO
238                  ENDDO
239               ELSE
240                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
241               ENDIF
242            END IF
243
244            IF ( ln_isf ) THEN   ! if ice shelf melting
245               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
246                  DO jj = 1, jpj
247                     DO ji = 1, jpi
248                        IF( misfkt(ji,jj) <= jk .AND. jk <= misfkb(ji,jj) ) THEN
249                           e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
250                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
251                        ENDIF
252                     END DO
253                  END DO
254               END DO
255            END IF
256            !
257            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
258               ! Before filtered scale factor at (u/v)-points
259               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
260               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
261               DO jk = 1, jpkm1
262                  DO jj = 1, jpj
263                     DO ji = 1, jpi
264                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
265                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
266                        !
267                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
268                        vb(ji,jj,jk) = zvf
269                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
270                        vn(ji,jj,jk) = va(ji,jj,jk)
271                     END DO
272                  END DO
273               END DO
274               !
275            ELSE                          ! Asselin filter applied on thickness weighted velocity
276               !
277               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
278               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
279               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
280               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
281               DO jk = 1, jpkm1
282                  DO jj = 1, jpj
283                     DO ji = 1, jpi                 
284                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
285                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
286                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
287                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
288                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
289                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
290                        !
291                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
292                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
293                        !
294                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
295                        vb(ji,jj,jk) = zvf
296                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
297                        vn(ji,jj,jk) = va(ji,jj,jk)
298                     END DO
299                  END DO
300               END DO
301               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
302               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
303               !
304               DEALLOCATE( ze3u_f , ze3v_f )
305            ENDIF
306            !
307         ENDIF
308         !
309         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
310            ! Revert "before" velocities to time split estimate
311            ! Doing it here also means that asselin filter contribution is removed 
312            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
313            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
314            DO jk = 2, jpkm1
315               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
316               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
317            END DO
318            DO jk = 1, jpkm1
319               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
320               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
321            END DO
322         ENDIF
323         !
324      ENDIF ! neuler =/0
325      !
326      ! Set "now" and "before" barotropic velocities for next time step:
327      ! JC: Would be more clever to swap variables than to make a full vertical
328      ! integration
329      !
330      !
331      IF(.NOT.ln_linssh ) THEN
332         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
333         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
334         DO jk = 2, jpkm1
335            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
336            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
337         END DO
338         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
339         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
340      ENDIF
341      !
342      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
343      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
344      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
345      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
346      DO jk = 2, jpkm1
347         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
348         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
349         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
350         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
351      END DO
352      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
353      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
354      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
355      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
356      !
357      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
358         CALL iom_put(  "ubar", un_b(:,:) )
359         CALL iom_put(  "vbar", vn_b(:,:) )
360      ENDIF
361      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
362         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
363         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
364         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
365      ENDIF
366      !
367      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
368         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
369      !
370      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
371      IF( l_trddyn     )   DEALLOCATE( zua, zva )
372      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
373      !
374   END SUBROUTINE dyn_nxt
375
376   !!=========================================================================
377END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.