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 NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90 @ 11027

Last change on this file since 11027 was 11027, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Final renaming conversions and removal of temporary pointers. All non-AGRIF SETTE tests are passing (including test cases). AGRIF tests compile and link but segment on first call to Agrif_Regrid. NST changes are therefore a work in progress but nothing is broken that was not broken before
  • Property svn:keywords set to Id
File size: 20.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 sbcrnf         ! river runoffs
31   USE sbcisf         ! ice shelf
32   USE phycst         ! physical constants
33   USE dynadv         ! dynamics: vector invariant versus flux form
34   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
35   USE domvvl         ! variable volume
36   USE bdy_oce   , ONLY: ln_bdy
37   USE bdydta         ! ocean open boundary conditions
38   USE bdydyn         ! ocean open boundary conditions
39   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
40   USE trd_oce        ! trends: ocean variables
41   USE trddyn         ! trend manager: dynamics
42   USE trdken         ! trend manager: kinetic energy
43   !
44   USE in_out_manager ! I/O manager
45   USE iom            ! I/O manager library
46   USE lbclnk         ! lateral boundary condition (or mpp link)
47   USE lib_mpp        ! MPP library
48   USE prtctl         ! Print control
49   USE timing         ! Timing
50#if defined key_agrif
51   USE agrif_oce_interp
52#endif   
53
54   IMPLICIT NONE
55   PRIVATE
56
57   PUBLIC    dyn_nxt   ! routine called by step.F90
58
59   !!----------------------------------------------------------------------
60   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
61   !! $Id$
62   !! Software governed by the CeCILL license (see ./LICENSE)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE dyn_nxt ( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa )
67      !!----------------------------------------------------------------------
68      !!                  ***  ROUTINE dyn_nxt  ***
69      !!                   
70      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
71      !!             condition on the after velocity, achieve the time stepping
72      !!             by applying the Asselin filter on now fields and swapping
73      !!             the fields.
74      !!
75      !! ** Method  : * Ensure after velocities transport matches time splitting
76      !!             estimate (ln_dynspg_ts=T)
77      !!
78      !!              * Apply lateral boundary conditions on after velocity
79      !!             at the local domain boundaries through lbc_lnk call,
80      !!             at the one-way open boundaries (ln_bdy=T),
81      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
82      !!
83      !!              * Apply the time filter applied and swap of the dynamics
84      !!             arrays to start the next time step:
85      !!                (puu(Kbb),pvv(Kbb)) = (puu(Kmm),pvv(Kmm)) + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Krhs),pvv(Krhs)) - 2 (puu(Kmm),pvv(Kmm)) ]
86      !!                (puu(Kmm),pvv(Kmm)) = (puu(Krhs),pvv(Krhs)).
87      !!             Note that with flux form advection and non linear free surface,
88      !!             the time filter is applied on thickness weighted velocity.
89      !!             As a result, dyn_nxt MUST be called after tra_nxt.
90      !!
91      !! ** Action :   puu(Kbb),pvv(Kbb)   filtered before horizontal velocity of next time-step
92      !!               puu(Kmm),pvv(Kmm)   now horizontal velocity of next time-step
93      !!----------------------------------------------------------------------
94      INTEGER, INTENT( in ) ::   kt                   ! ocean time-step index
95      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs, Kaa  ! time level indices
96      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv        !  velocity arrays
97      !
98      INTEGER  ::   ji, jj, jk   ! dummy loop indices
99      INTEGER  ::   ikt          ! local integers
100      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
101      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
102      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
103      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
104      !!----------------------------------------------------------------------
105      !
106      IF( ln_timing    )   CALL timing_start('dyn_nxt')
107      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
108      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
109      !
110      IF( kt == nit000 ) THEN
111         IF(lwp) WRITE(numout,*)
112         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
113         IF(lwp) WRITE(numout,*) '~~~~~~~'
114      ENDIF
115
116      IF ( ln_dynspg_ts ) THEN
117         ! Ensure below that barotropic velocities match time splitting estimate
118         ! Compute actual transport and replace it with ts estimate at "after" time step
119         zue(:,:) = e3u(:,:,1,Krhs) * puu(:,:,1,Krhs) * umask(:,:,1)
120         zve(:,:) = e3v(:,:,1,Krhs) * pvv(:,:,1,Krhs) * vmask(:,:,1)
121         DO jk = 2, jpkm1
122            zue(:,:) = zue(:,:) + e3u(:,:,jk,Krhs) * puu(:,:,jk,Krhs) * umask(:,:,jk)
123            zve(:,:) = zve(:,:) + e3v(:,:,jk,Krhs) * pvv(:,:,jk,Krhs) * vmask(:,:,jk)
124         END DO
125         DO jk = 1, jpkm1
126            puu(:,:,jk,Krhs) = ( puu(:,:,jk,Krhs) - zue(:,:) * r1_hu_a(:,:) + uu_b(:,:,Krhs) ) * umask(:,:,jk)
127            pvv(:,:,jk,Krhs) = ( pvv(:,:,jk,Krhs) - zve(:,:) * r1_hv_a(:,:) + vv_b(:,:,Krhs) ) * vmask(:,:,jk)
128         END DO
129         !
130         IF( .NOT.ln_bt_fw ) THEN
131            ! Remove advective velocity from "now velocities"
132            ! prior to asselin filtering     
133            ! In the forward case, this is done below after asselin filtering   
134            ! so that asselin contribution is removed at the same time
135            DO jk = 1, jpkm1
136               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu_n(:,:) + uu_b(:,:,Kmm) )*umask(:,:,jk)
137               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv_n(:,:) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
138            END DO 
139         ENDIF
140      ENDIF
141
142      ! Update after velocity on domain lateral boundaries
143      ! --------------------------------------------------     
144# if defined key_agrif
145      Krhs_a = Krhs ; CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
146# endif
147      !
148      CALL lbc_lnk_multi( 'dynnxt', puu(:,:,:,Krhs), 'U', -1., pvv(:,:,:,Krhs), 'V', -1. )     !* local domain boundaries
149      !
150      !                                !* BDY open boundaries
151      !! IMMERSE development : Following the general pattern for the new code we want to pass in the
152      !!                       velocities to bdy_dyn as arguments so here we use "uu" and "vv" even
153      !!                       though we haven't converted the velocity names in the rest of dynnxt.F90
154      !!                       because it will be completely rewritten. DS.
155      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, uu, vv, Kaa )
156      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, uu, vv, Kaa, dyn3d_only=.true. )
157
158!!$   Do we need a call to bdy_vol here??
159      !
160      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
161         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
162         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
163         !
164         !                                  ! Kinetic energy and Conversion
165         IF( ln_KE_trd  )   CALL trd_dyn( puu(:,:,:,Krhs), pvv(:,:,:,Krhs), jpdyn_ken, kt, Kmm )
166         !
167         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
168            zua(:,:,:) = ( puu(:,:,:,Krhs) - puu(:,:,:,Kbb) ) * z1_2dt
169            zva(:,:,:) = ( pvv(:,:,:,Krhs) - pvv(:,:,:,Kbb) ) * z1_2dt
170            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
171            CALL iom_put( "vtrd_tot", zva )
172         ENDIF
173         !
174         zua(:,:,:) = puu(:,:,:,Kmm)             ! save the now velocity before the asselin filter
175         zva(:,:,:) = pvv(:,:,:,Kmm)             ! (caution: there will be a shift by 1 timestep in the
176         !                                  !  computation of the asselin filter trends)
177      ENDIF
178
179      ! Time filter and swap of dynamics arrays
180      ! ------------------------------------------
181      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
182         DO jk = 1, jpkm1
183            puu(:,:,jk,Kmm) = puu(:,:,jk,Krhs)                         ! puu(:,:,:,Kmm) <-- puu(:,:,:,Krhs)
184            pvv(:,:,jk,Kmm) = pvv(:,:,jk,Krhs)
185         END DO
186         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
187!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a 
188            DO jk = 1, jpkm1
189!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
190!               e3u(:,:,jk,Kbb) = e3u(:,:,jk,Kmm)
191!               e3v(:,:,jk,Kbb) = e3v(:,:,jk,Kmm)
192               !
193               e3t(:,:,jk,Kmm) = e3t(:,:,jk,Krhs)
194               e3u(:,:,jk,Kmm) = e3u(:,:,jk,Krhs)
195               e3v(:,:,jk,Kmm) = e3v(:,:,jk,Krhs)
196            END DO
197!!gm BUG end
198         ENDIF
199                                                            !
200         
201      ELSE                                             !* Leap-Frog : Asselin filter and swap
202         !                                ! =============!
203         IF( ln_linssh ) THEN             ! Fixed volume !
204            !                             ! =============!
205            DO jk = 1, jpkm1                             
206               DO jj = 1, jpj
207                  DO ji = 1, jpi   
208                     zuf = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Krhs) )
209                     zvf = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Krhs) )
210                     !
211                     puu(ji,jj,jk,Kbb) = zuf                      ! puu(:,:,:,Kbb) <-- filtered velocity
212                     pvv(ji,jj,jk,Kbb) = zvf
213                     puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Krhs)             ! puu(:,:,:,Kmm) <-- puu(:,:,:,Krhs)
214                     pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Krhs)
215                  END DO
216               END DO
217            END DO
218            !                             ! ================!
219         ELSE                             ! Variable volume !
220            !                             ! ================!
221            ! Before scale factor at t-points
222            ! (used as a now filtered scale factor until the swap)
223            ! ----------------------------------------------------
224            DO jk = 1, jpkm1
225               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) + atfp * ( e3t(:,:,jk,Kbb) - 2._wp * e3t(:,:,jk,Kmm) + e3t(:,:,jk,Krhs) )
226            END DO
227            ! Add volume filter correction: compatibility with tracer advection scheme
228            ! => time filter + conservation correction (only at the first level)
229            zcoef = atfp * rdt * r1_rau0
230
231            e3t(:,:,1,Kbb) = e3t(:,:,1,Kbb) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
232
233            IF ( ln_rnf ) THEN
234               IF( ln_rnf_depth ) THEN
235                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
236                     DO jj = 1, jpj
237                        DO ji = 1, jpi
238                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
239                               e3t(ji,jj,jk,Kbb) =   e3t(ji,jj,jk,Kbb) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
240                                      &          * ( e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
241                           ENDIF
242                        ENDDO
243                     ENDDO
244                  ENDDO
245               ELSE
246                  e3t(:,:,1,Kbb) = e3t(:,:,1,Kbb) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
247               ENDIF
248            END IF
249
250            IF ( ln_isf ) THEN   ! if ice shelf melting
251               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
252                  DO jj = 1, jpj
253                     DO ji = 1, jpi
254                        IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN
255                           e3t(ji,jj,jk,Kbb) = e3t(ji,jj,jk,Kbb) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
256                                &          * ( e3t(ji,jj,jk,Kmm) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
257                        ELSEIF ( jk==misfkb(ji,jj) ) THEN
258                           e3t(ji,jj,jk,Kbb) = e3t(ji,jj,jk,Kbb) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
259                                &          * ( e3t(ji,jj,jk,Kmm) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
260                        ENDIF
261                     END DO
262                  END DO
263               END DO
264            END IF
265            !
266            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
267               ! Before filtered scale factor at (u/v)-points
268               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
269               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
270               DO jk = 1, jpkm1
271                  DO jj = 1, jpj
272                     DO ji = 1, jpi
273                        zuf = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Krhs) )
274                        zvf = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Krhs) )
275                        !
276                        puu(ji,jj,jk,Kbb) = zuf                      ! puu(:,:,:,Kbb) <-- filtered velocity
277                        pvv(ji,jj,jk,Kbb) = zvf
278                        puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Krhs)             ! puu(:,:,:,Kmm) <-- puu(:,:,:,Krhs)
279                        pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Krhs)
280                     END DO
281                  END DO
282               END DO
283               !
284            ELSE                          ! Asselin filter applied on thickness weighted velocity
285               !
286               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
287               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
288               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), ze3u_f, 'U' )
289               CALL dom_vvl_interpol( e3t(:,:,:,Kbb), ze3v_f, 'V' )
290               DO jk = 1, jpkm1
291                  DO jj = 1, jpj
292                     DO ji = 1, jpi                 
293                        zue3a = e3u(ji,jj,jk,Krhs) * puu(ji,jj,jk,Krhs)
294                        zve3a = e3v(ji,jj,jk,Krhs) * pvv(ji,jj,jk,Krhs)
295                        zue3n = e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
296                        zve3n = e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
297                        zue3b = e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb)
298                        zve3b = e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb)
299                        !
300                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
301                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
302                        !
303                        puu(ji,jj,jk,Kbb) = zuf                     ! puu(:,:,:,Kbb) <-- filtered velocity
304                        pvv(ji,jj,jk,Kbb) = zvf
305                        puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Krhs)            ! puu(:,:,:,Kmm) <-- puu(:,:,:,Krhs)
306                        pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Krhs)
307                     END DO
308                  END DO
309               END DO
310               e3u(:,:,1:jpkm1,Kbb) = ze3u_f(:,:,1:jpkm1)        ! e3u(:,:,:,Kbb) <-- filtered scale factor
311               e3v(:,:,1:jpkm1,Kbb) = ze3v_f(:,:,1:jpkm1)
312               !
313               DEALLOCATE( ze3u_f , ze3v_f )
314            ENDIF
315            !
316         ENDIF
317         !
318         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
319            ! Revert "before" velocities to time split estimate
320            ! Doing it here also means that asselin filter contribution is removed 
321            zue(:,:) = e3u(:,:,1,Kbb) * puu(:,:,1,Kbb) * umask(:,:,1)
322            zve(:,:) = e3v(:,:,1,Kbb) * pvv(:,:,1,Kbb) * vmask(:,:,1)   
323            DO jk = 2, jpkm1
324               zue(:,:) = zue(:,:) + e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb) * umask(:,:,jk)
325               zve(:,:) = zve(:,:) + e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb) * vmask(:,:,jk)   
326            END DO
327            DO jk = 1, jpkm1
328               puu(:,:,jk,Kbb) = puu(:,:,jk,Kbb) - (zue(:,:) * r1_hu_n(:,:) - uu_b(:,:,Kmm)) * umask(:,:,jk)
329               pvv(:,:,jk,Kbb) = pvv(:,:,jk,Kbb) - (zve(:,:) * r1_hv_n(:,:) - vv_b(:,:,Kmm)) * vmask(:,:,jk)
330            END DO
331         ENDIF
332         !
333      ENDIF ! neuler =/0
334      !
335      ! Set "now" and "before" barotropic velocities for next time step:
336      ! JC: Would be more clever to swap variables than to make a full vertical
337      ! integration
338      !
339      !
340      IF(.NOT.ln_linssh ) THEN
341         hu_b(:,:) = e3u(:,:,1,Kbb) * umask(:,:,1)
342         hv_b(:,:) = e3v(:,:,1,Kbb) * vmask(:,:,1)
343         DO jk = 2, jpkm1
344            hu_b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
345            hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
346         END DO
347         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
348         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
349      ENDIF
350      !
351      uu_b(:,:,Kmm) = e3u(:,:,1,Krhs) * puu(:,:,1,Kmm) * umask(:,:,1)
352      uu_b(:,:,Kbb) = e3u(:,:,1,Kbb) * puu(:,:,1,Kbb) * umask(:,:,1)
353      vv_b(:,:,Kmm) = e3v(:,:,1,Krhs) * pvv(:,:,1,Kmm) * vmask(:,:,1)
354      vv_b(:,:,Kbb) = e3v(:,:,1,Kbb) * pvv(:,:,1,Kbb) * vmask(:,:,1)
355      DO jk = 2, jpkm1
356         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + e3u(:,:,jk,Krhs) * puu(:,:,jk,Kmm) * umask(:,:,jk)
357         uu_b(:,:,Kbb) = uu_b(:,:,Kbb) + e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb) * umask(:,:,jk)
358         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + e3v(:,:,jk,Krhs) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
359         vv_b(:,:,Kbb) = vv_b(:,:,Kbb) + e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb) * vmask(:,:,jk)
360      END DO
361      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu_a(:,:)
362      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv_a(:,:)
363      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu_b(:,:)
364      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv_b(:,:)
365      !
366      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
367         CALL iom_put(  "ubar", uu_b(:,:,Kmm) )
368         CALL iom_put(  "vbar", vv_b(:,:,Kmm) )
369      ENDIF
370      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
371         zua(:,:,:) = ( puu(:,:,:,Kbb) - zua(:,:,:) ) * z1_2dt
372         zva(:,:,:) = ( pvv(:,:,:,Kbb) - zva(:,:,:) ) * z1_2dt
373         CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm )
374      ENDIF
375      !
376      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kmm), clinfo1=' nxt  - Un: ', mask1=umask,   &
377         &                       tab3d_2=pvv(:,:,:,Kmm), clinfo2=' Vn: '       , mask2=vmask )
378      !
379      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
380      IF( l_trddyn     )   DEALLOCATE( zua, zva )
381      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
382      !
383   END SUBROUTINE dyn_nxt
384
385   !!=========================================================================
386END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.