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 @ 7513

Last change on this file since 7513 was 7513, checked in by mocavero, 7 years ago

delete optional ompdo statements
  • Property svn:keywords set to Id
File size: 21.4 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
118!$OMP DO schedule(static) private(jj, ji)
119         DO jj = 1, jpj
120            DO ji = 1, jpi
121               zue(ji,jj) = e3u_a(ji,jj,1) * ua(ji,jj,1) * umask(ji,jj,1)
122               zve(ji,jj) = e3v_a(ji,jj,1) * va(ji,jj,1) * vmask(ji,jj,1)
123            END DO
124         END DO
125         DO jk = 2, jpkm1
126!$OMP DO schedule(static) private(jj,ji)
127            DO jj = 1, jpj
128               DO ji = 1, jpi
129                  zue(ji,jj) = zue(ji,jj) + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
130                  zve(ji,jj) = zve(ji,jj) + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
131               END DO
132            END DO
133         END DO
134!$OMP DO schedule(static) private(jk)
135         DO jk = 1, jpkm1
136            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
137            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
138         END DO
139!$OMP END DO NOWAIT
140!$OMP END PARALLEL
141         !
142         IF( .NOT.ln_bt_fw ) THEN
143            ! Remove advective velocity from "now velocities"
144            ! prior to asselin filtering     
145            ! In the forward case, this is done below after asselin filtering   
146            ! so that asselin contribution is removed at the same time
147!$OMP PARALLEL DO schedule(static) private(jk)
148            DO jk = 1, jpkm1
149               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
150               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
151            END DO 
152         ENDIF
153      ENDIF
154
155      ! Update after velocity on domain lateral boundaries
156      ! --------------------------------------------------     
157# if defined key_agrif
158      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
159# endif
160      !
161      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
162      CALL lbc_lnk( va, 'V', -1. ) 
163      !
164# if defined key_bdy
165      !                                !* BDY open boundaries
166      IF( lk_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
167      IF( lk_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
168
169!!$   Do we need a call to bdy_vol here??
170      !
171# endif
172      !
173      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
174         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
175         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
176         !
177         !                                  ! Kinetic energy and Conversion
178         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
179         !
180         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
181!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
182         DO jk = 1, jpk
183            DO jj = 1, jpj
184               DO ji = 1, jpi
185                  zua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_2dt
186                  zva(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_2dt
187               END DO
188            END DO
189         END DO
190         CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
191         CALL iom_put( "vtrd_tot", zva )
192         ENDIF
193         !
194!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
195         DO jk = 1, jpk
196            DO jj = 1, jpj
197               DO ji = 1, jpi
198                  zua(ji,jj,jk) = un(ji,jj,jk)             ! save the now velocity before the asselin filter
199                  zva(ji,jj,jk) = vn(ji,jj,jk)             ! (caution: there will be a shift by 1 timestep in the
200               END DO
201            END DO
202         END DO
203         !                                  !  computation of the asselin filter trends)
204      ENDIF
205
206      ! Time filter and swap of dynamics arrays
207      ! ------------------------------------------
208      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
209!$OMP PARALLEL
210!$OMP DO schedule(static) private(jk)
211         DO jk = 1, jpkm1
212            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
213            vn(:,:,jk) = va(:,:,jk)
214         END DO
215!$OMP END DO NOWAIT
216         IF(.NOT.ln_linssh ) THEN
217!$OMP DO schedule(static) private(jk)
218            DO jk = 1, jpkm1
219               e3t_b(:,:,jk) = e3t_n(:,:,jk)
220               e3u_b(:,:,jk) = e3u_n(:,:,jk)
221               e3v_b(:,:,jk) = e3v_n(:,:,jk)
222            END DO
223!$OMP END DO NOWAIT
224         ENDIF
225!$OMP END PARALLEL
226      ELSE                                             !* Leap-Frog : Asselin filter and swap
227         !                                ! =============!
228         IF( ln_linssh ) THEN             ! Fixed volume !
229            !                             ! =============!
230!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zuf, zvf)
231            DO jk = 1, jpkm1                             
232               DO jj = 1, jpj
233                  DO ji = 1, jpi   
234                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
235                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
236                     !
237                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
238                     vb(ji,jj,jk) = zvf
239                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
240                     vn(ji,jj,jk) = va(ji,jj,jk)
241                  END DO
242               END DO
243            END DO
244            !                             ! ================!
245         ELSE                             ! Variable volume !
246            !                             ! ================!
247            ! Before scale factor at t-points
248            ! (used as a now filtered scale factor until the swap)
249            ! ----------------------------------------------------
250            IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN    ! No asselin filtering on thicknesses if forward time splitting
251!$OMP PARALLEL DO schedule(static) private(jj,ji)
252               DO jj = 1, jpj
253                  DO ji = 1, jpi
254                     e3t_b(ji,jj,1:jpkm1) = e3t_n(ji,jj,1:jpkm1)
255                  END DO
256               END DO
257            ELSE
258!$OMP PARALLEL DO schedule(static) private(jk)
259               DO jk = 1, jpkm1
260                  e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
261               END DO
262               ! Add volume filter correction: compatibility with tracer advection scheme
263               ! => time filter + conservation correction (only at the first level)
264               zcoef = atfp * rdt * r1_rau0
265               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting
266!$OMP PARALLEL DO schedule(static) private(jj,ji)
267                  DO jj = 1, jpj
268                     DO ji = 1, jpi
269                        e3t_b(ji,jj,1) = e3t_b(ji,jj,1) - zcoef * ( emp_b(ji,jj) - emp(ji,jj) &
270                                 &                      - rnf_b(ji,jj) + rnf(ji,jj) ) * tmask(ji,jj,1)
271                     END DO
272                  END DO
273               ELSE                     ! if ice shelf melting
274!$OMP PARALLEL DO schedule(static) private(jj,ji,ikt)
275                  DO jj = 1, jpj
276                     DO ji = 1, jpi
277                        ikt = mikt(ji,jj)
278                        e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  &
279                           &                                           - rnf_b   (ji,jj) + rnf   (ji,jj)  &
280                           &                                           + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,ikt)
281                     END DO
282                  END DO
283               END IF
284            ENDIF
285            !
286            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
287               ! Before filtered scale factor at (u/v)-points
288               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
289               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
290!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zuf, zvf)
291               DO jk = 1, jpkm1
292                  DO jj = 1, jpj
293                     DO ji = 1, jpi
294                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
295                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
296                        !
297                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
298                        vb(ji,jj,jk) = zvf
299                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
300                        vn(ji,jj,jk) = va(ji,jj,jk)
301                     END DO
302                  END DO
303               END DO
304               !
305            ELSE                          ! Asselin filter applied on thickness weighted velocity
306               !
307               CALL wrk_alloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
308               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
309               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
310               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
311!$OMP PARALLEL
312!$OMP DO schedule(static) private(jk, jj, ji, zue3a, zve3a, zue3n, zve3n, zue3b, zve3b, zuf, zvf)
313               DO jk = 1, jpkm1
314                  DO jj = 1, jpj
315                     DO ji = 1, jpi                 
316                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
317                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
318                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
319                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
320                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
321                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
322                        !
323                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
324                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
325                        !
326                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
327                        vb(ji,jj,jk) = zvf
328                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
329                        vn(ji,jj,jk) = va(ji,jj,jk)
330                     END DO
331                  END DO
332               END DO
333!$OMP DO schedule(static) private(jj, ji)
334                  DO jj = 1, jpj
335                     DO ji = 1, jpi
336                        e3u_b(ji,jj,1:jpkm1) = ze3u_f(ji,jj,1:jpkm1)        ! e3u_b <-- filtered scale factor
337                        e3v_b(ji,jj,1:jpkm1) = ze3v_f(ji,jj,1:jpkm1)
338                     END DO
339                  END DO
340!$OMP END PARALLEL
341               !
342               CALL wrk_dealloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
343            ENDIF
344            !
345         ENDIF
346         !
347         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
348            ! Revert "before" velocities to time split estimate
349            ! Doing it here also means that asselin filter contribution is removed 
350!$OMP PARALLEL
351!$OMP DO schedule(static) private(jj, ji)
352            DO jj = 1, jpj
353               DO ji = 1, jpi
354                  zue(ji,jj) = e3u_b(ji,jj,1) * ub(ji,jj,1) * umask(ji,jj,1)
355                  zve(ji,jj) = e3v_b(ji,jj,1) * vb(ji,jj,1) * vmask(ji,jj,1)   
356               END DO
357            END DO
358            DO jk = 2, jpkm1
359!$OMP DO schedule(static) private(jj, ji)
360               DO jj = 1, jpj
361                  DO ji = 1, jpi
362                     zue(ji,jj) = zue(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
363                     zve(ji,jj) = zve(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)   
364                  END DO
365               END DO
366            END DO
367!$OMP DO schedule(static) private(jk)
368            DO jk = 1, jpkm1
369               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
370               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
371            END DO
372!$OMP END DO NOWAIT
373!$OMP END PARALLEL
374         ENDIF
375         !
376      ENDIF ! neuler =/0
377      !
378      ! Set "now" and "before" barotropic velocities for next time step:
379      ! JC: Would be more clever to swap variables than to make a full vertical
380      ! integration
381      !
382      !
383      IF(.NOT.ln_linssh ) THEN
384!$OMP PARALLEL
385!$OMP DO schedule(static) private(jj, ji)
386         DO jj = 1, jpj
387            DO ji = 1, jpi
388               hu_b(ji,jj) = e3u_b(ji,jj,1) * umask(ji,jj,1)
389               hv_b(ji,jj) = e3v_b(ji,jj,1) * vmask(ji,jj,1)
390            END DO
391         END DO
392         DO jk = 2, jpkm1
393!$OMP DO schedule(static) private(jj, ji)
394            DO jj = 1, jpj
395               DO ji = 1, jpi
396                  hu_b(ji,jj) = hu_b(ji,jj) + e3u_b(ji,jj,jk) * umask(ji,jj,jk)
397                  hv_b(ji,jj) = hv_b(ji,jj) + e3v_b(ji,jj,jk) * vmask(ji,jj,jk)
398               END DO
399            END DO
400         END DO
401!$OMP DO schedule(static) private(jj, ji)
402         DO jj = 1, jpj
403            DO ji = 1, jpi
404               r1_hu_b(ji,jj) = ssumask(ji,jj) / ( hu_b(ji,jj) + 1._wp - ssumask(ji,jj) )
405               r1_hv_b(ji,jj) = ssvmask(ji,jj) / ( hv_b(ji,jj) + 1._wp - ssvmask(ji,jj) )
406            END DO
407         END DO
408!$OMP END PARALLEL
409      ENDIF
410      !
411!$OMP PARALLEL
412!$OMP DO schedule(static) private(jj, ji)
413      DO jj = 1, jpj
414         DO ji = 1, jpi
415            un_b(ji,jj) = e3u_a(ji,jj,1) * un(ji,jj,1) * umask(ji,jj,1)
416            ub_b(ji,jj) = e3u_b(ji,jj,1) * ub(ji,jj,1) * umask(ji,jj,1)
417            vn_b(ji,jj) = e3v_a(ji,jj,1) * vn(ji,jj,1) * vmask(ji,jj,1)
418            vb_b(ji,jj) = e3v_b(ji,jj,1) * vb(ji,jj,1) * vmask(ji,jj,1)
419         END DO
420      END DO
421      DO jk = 2, jpkm1
422!$OMP DO schedule(static) private(jj, ji)
423         DO jj = 1, jpj
424            DO ji = 1, jpi
425               un_b(ji,jj) = un_b(ji,jj) + e3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
426               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
427               vn_b(ji,jj) = vn_b(ji,jj) + e3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
428               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
429            END DO
430         END DO
431      END DO
432!$OMP DO schedule(static) private(jj, ji)
433      DO jj = 1, jpj
434         DO ji = 1, jpi
435            un_b(ji,jj) = un_b(ji,jj) * r1_hu_a(ji,jj)
436            vn_b(ji,jj) = vn_b(ji,jj) * r1_hv_a(ji,jj)
437            ub_b(ji,jj) = ub_b(ji,jj) * r1_hu_b(ji,jj)
438            vb_b(ji,jj) = vb_b(ji,jj) * r1_hv_b(ji,jj)
439         END DO
440      END DO
441!$OMP END PARALLEL
442      !
443      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
444         CALL iom_put(  "ubar", un_b(:,:) )
445         CALL iom_put(  "vbar", vn_b(:,:) )
446      ENDIF
447      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
448!$OMP DO schedule(static) private(jk, jj, ji)
449         DO jk = 1, jpkm1
450            DO jj = 1, jpj
451               DO ji = 1, jpi
452                  zua(ji,jj,jk) = ( ub(ji,jj,jk) - zua(ji,jj,jk) ) * z1_2dt
453                  zva(ji,jj,jk) = ( vb(ji,jj,jk) - zva(ji,jj,jk) ) * z1_2dt
454               END DO
455            END DO
456        END DO
457        CALL trd_dyn( zua, zva, jpdyn_atf, kt )
458      ENDIF
459      !
460      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
461         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
462      !
463      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj,       zue, zve )
464      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk,   zua, zva )
465      !
466      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
467      !
468   END SUBROUTINE dyn_nxt
469
470   !!=========================================================================
471END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.