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

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

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