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/OPA_SRC/DYN – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 9090

Last change on this file since 9090 was 9090, checked in by flavoni, 6 years ago

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