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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 7698

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

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