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

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 6748

Last change on this file since 6748 was 6748, checked in by mocavero, 8 years ago

GYRE hybrid parallelization
  • Property svn:keywords set to Id
File size: 18.6 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 phycst         ! physical constants
31   USE dynadv         ! dynamics: vector invariant versus flux form
32   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
33   USE domvvl         ! variable volume
34   USE bdy_oce        ! ocean open boundary conditions
35   USE bdydta         ! ocean open boundary conditions
36   USE bdydyn         ! ocean open boundary conditions
37   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
38   USE trd_oce        ! trends: ocean variables
39   USE trddyn         ! trend manager: dynamics
40   USE trdken         ! trend manager: kinetic energy
41   !
42   USE in_out_manager ! I/O manager
43   USE iom            ! I/O manager library
44   USE lbclnk         ! lateral boundary condition (or mpp link)
45   USE lib_mpp        ! MPP library
46   USE wrk_nemo       ! Memory Allocation
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 3.3 , NEMO Consortium (2010)
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 (lk_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), POINTER, DIMENSION(:,:)   ::  zue, zve
100      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
101      !!----------------------------------------------------------------------
102      !
103      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
104      !
105      IF( ln_dynspg_ts       )   CALL wrk_alloc( jpi,jpj,       zue, zve)
106      IF( l_trddyn           )   CALL wrk_alloc( jpi,jpj,jpk,   zua, zva)
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!$OMP PARALLEL WORKSHARE
118         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
119         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
120!$OMP END PARALLEL WORKSHARE
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!$OMP PARALLEL DO schedule(static) private(jk)
126         DO jk = 1, jpkm1
127            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
128            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
129         END DO
130         !
131         IF( .NOT.ln_bt_fw ) THEN
132            ! Remove advective velocity from "now velocities"
133            ! prior to asselin filtering     
134            ! In the forward case, this is done below after asselin filtering   
135            ! so that asselin contribution is removed at the same time
136!$OMP PARALLEL DO schedule(static) private(jk)
137            DO jk = 1, jpkm1
138               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
139               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
140            END DO 
141         ENDIF
142      ENDIF
143
144      ! Update after velocity on domain lateral boundaries
145      ! --------------------------------------------------     
146# if defined key_agrif
147      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
148# endif
149      !
150      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
151      CALL lbc_lnk( va, 'V', -1. ) 
152      !
153# if defined key_bdy
154      !                                !* BDY open boundaries
155      IF( lk_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
156      IF( lk_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
157
158!!$   Do we need a call to bdy_vol here??
159      !
160# endif
161      !
162      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
163         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
164         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
165         !
166         !                                  ! Kinetic energy and Conversion
167         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
168         !
169         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
170!$OMP PARALLEL WORKSHARE
171            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
172            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
173!$OMP END PARALLEL WORKSHARE
174            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
175            CALL iom_put( "vtrd_tot", zva )
176         ENDIF
177         !
178!$OMP PARALLEL WORKSHARE
179         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
180         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
181!$OMP END PARALLEL WORKSHARE
182         !                                  !  computation of the asselin filter trends)
183      ENDIF
184
185      ! Time filter and swap of dynamics arrays
186      ! ------------------------------------------
187      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
188!$OMP PARALLEL DO schedule(static) private(jk)
189         DO jk = 1, jpkm1
190            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
191            vn(:,:,jk) = va(:,:,jk)
192         END DO
193         IF(.NOT.ln_linssh ) THEN
194!$OMP PARALLEL DO schedule(static) private(jk)
195            DO jk = 1, jpkm1
196               e3t_b(:,:,jk) = e3t_n(:,:,jk)
197               e3u_b(:,:,jk) = e3u_n(:,:,jk)
198               e3v_b(:,:,jk) = e3v_n(:,:,jk)
199            END DO
200         ENDIF
201      ELSE                                             !* Leap-Frog : Asselin filter and swap
202         !                                ! =============!
203         IF( ln_linssh ) THEN             ! Fixed volume !
204            !                             ! =============!
205!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zuf, zvf)
206            DO jk = 1, jpkm1                             
207               DO jj = 1, jpj
208                  DO ji = 1, jpi   
209                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
210                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
211                     !
212                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
213                     vb(ji,jj,jk) = zvf
214                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
215                     vn(ji,jj,jk) = va(ji,jj,jk)
216                  END DO
217               END DO
218            END DO
219            !                             ! ================!
220         ELSE                             ! Variable volume !
221            !                             ! ================!
222            ! Before scale factor at t-points
223            ! (used as a now filtered scale factor until the swap)
224            ! ----------------------------------------------------
225            IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN    ! No asselin filtering on thicknesses if forward time splitting
226               e3t_b(:,:,1:jpkm1) = e3t_n(:,:,1:jpkm1)
227            ELSE
228!$OMP PARALLEL DO schedule(static) private(jk)
229               DO jk = 1, jpkm1
230                  e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
231               END DO
232               ! Add volume filter correction: compatibility with tracer advection scheme
233               ! => time filter + conservation correction (only at the first level)
234               zcoef = atfp * rdt * r1_rau0
235               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting
236                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) &
237                                 &                      - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
238               ELSE                     ! if ice shelf melting
239                  DO jj = 1, jpj
240                     DO ji = 1, jpi
241                        ikt = mikt(ji,jj)
242                        e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  &
243                           &                                           - rnf_b   (ji,jj) + rnf   (ji,jj)  &
244                           &                                           + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,ikt)
245                     END DO
246                  END DO
247               END IF
248            ENDIF
249            !
250            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
251               ! Before filtered scale factor at (u/v)-points
252               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
253               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
254!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zuf, zvf)
255               DO jk = 1, jpkm1
256                  DO jj = 1, jpj
257                     DO ji = 1, jpi
258                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
259                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
260                        !
261                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
262                        vb(ji,jj,jk) = zvf
263                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
264                        vn(ji,jj,jk) = va(ji,jj,jk)
265                     END DO
266                  END DO
267               END DO
268               !
269            ELSE                          ! Asselin filter applied on thickness weighted velocity
270               !
271               CALL wrk_alloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
272               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
273               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
274               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
275!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zue3a, zve3a, zue3n, zve3n, zue3b, zve3b, zuf, zvf)
276               DO jk = 1, jpkm1
277                  DO jj = 1, jpj
278                     DO ji = 1, jpi                 
279                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
280                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
281                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
282                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
283                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
284                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
285                        !
286                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
287                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
288                        !
289                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
290                        vb(ji,jj,jk) = zvf
291                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
292                        vn(ji,jj,jk) = va(ji,jj,jk)
293                     END DO
294                  END DO
295               END DO
296               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
297               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
298               !
299               CALL wrk_dealloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
300            ENDIF
301            !
302         ENDIF
303         !
304         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
305            ! Revert "before" velocities to time split estimate
306            ! Doing it here also means that asselin filter contribution is removed 
307            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
308            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
309            DO jk = 2, jpkm1
310               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
311               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
312            END DO
313!$OMP PARALLEL DO schedule(static) private(jk)
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!$OMP PARALLEL WORKSHARE
339      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
340      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
341      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
342      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
343!$OMP END PARALLEL WORKSHARE
344      DO jk = 2, jpkm1
345         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
346         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
347         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
348         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
349      END DO
350!$OMP PARALLEL WORKSHARE
351      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
352      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
353      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
354      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
355!$OMP END PARALLEL WORKSHARE
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!$OMP PARALLEL WORKSHARE
363         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
364         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
365!$OMP END PARALLEL WORKSHARE
366         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
367      ENDIF
368      !
369      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
370         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
371      !
372      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj,       zue, zve )
373      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk,   zua, zva )
374      !
375      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
376      !
377   END SUBROUTINE dyn_nxt
378
379   !!=========================================================================
380END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.