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/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/DYN/dynnxt.F90 @ 13454

Last change on this file since 13454 was 13454, checked in by dancopsey, 4 years ago

Merge in Catherine Guiavarch's change to implement implicit ice drag to improve stability. Her changes were implemented in:

https://code.metoffice.gov.uk/trac/roses-u/changeset?reponame=&new=161415%40b%2Fv%2F1%2F1%2F4&old=161163%40b%2Fv%2F1%2F1%2F4

File size: 20.1 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   USE zdfdrg, ONLY: ln_drgice_imp, rCdU_top
51#if defined key_agrif
52   USE agrif_oce_interp
53#endif
54
55   IMPLICIT NONE
56   PRIVATE
57
58   PUBLIC    dyn_nxt   ! routine called by step.F90
59
60   !!----------------------------------------------------------------------
61   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
62   !! $Id$
63   !! Software governed by the CeCILL license (see ./LICENSE)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE dyn_nxt ( kt )
68      !!----------------------------------------------------------------------
69      !!                  ***  ROUTINE dyn_nxt  ***
70      !!                   
71      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
72      !!             condition on the after velocity, achieve the time stepping
73      !!             by applying the Asselin filter on now fields and swapping
74      !!             the fields.
75      !!
76      !! ** Method  : * Ensure after velocities transport matches time splitting
77      !!             estimate (ln_dynspg_ts=T)
78      !!
79      !!              * Apply lateral boundary conditions on after velocity
80      !!             at the local domain boundaries through lbc_lnk call,
81      !!             at the one-way open boundaries (ln_bdy=T),
82      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
83      !!
84      !!              * Apply the time filter applied and swap of the dynamics
85      !!             arrays to start the next time step:
86      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
87      !!                (un,vn) = (ua,va).
88      !!             Note that with flux form advection and non linear free surface,
89      !!             the time filter is applied on thickness weighted velocity.
90      !!             As a result, dyn_nxt MUST be called after tra_nxt.
91      !!
92      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
93      !!               un,vn   now horizontal velocity of next time-step
94      !!----------------------------------------------------------------------
95      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
96      !
97      INTEGER  ::   ji, jj, jk   ! dummy loop indices
98      INTEGER  ::   ikt          ! local integers
99      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
100      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
101      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
102      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zutau, zvtau
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_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
120         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
121         DO jk = 2, jpkm1
122            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
123            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
124         END DO
125         DO jk = 1, jpkm1
126            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
127            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * 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               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
137               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*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      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
146# endif
147      !
148      CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
149      !
150      !                                !* BDY open boundaries
151      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
152      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
153
154!!$   Do we need a call to bdy_vol here??
155      !
156      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
157         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
158         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
159         !
160         !                                  ! Kinetic energy and Conversion
161         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
162         !
163         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
164            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
165            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
166            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
167            CALL iom_put( "vtrd_tot", zva )
168         ENDIF
169         !
170         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
171         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
172         !                                  !  computation of the asselin filter trends)
173      ENDIF
174
175      ! Time filter and swap of dynamics arrays
176      ! ------------------------------------------
177      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
178         DO jk = 1, jpkm1
179            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
180            vn(:,:,jk) = va(:,:,jk)
181         END DO
182         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
183!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a 
184            DO jk = 1, jpkm1
185!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
186!               e3u_b(:,:,jk) = e3u_n(:,:,jk)
187!               e3v_b(:,:,jk) = e3v_n(:,:,jk)
188               !
189               e3t_n(:,:,jk) = e3t_a(:,:,jk)
190               e3u_n(:,:,jk) = e3u_a(:,:,jk)
191               e3v_n(:,:,jk) = e3v_a(:,:,jk)
192            END DO
193!!gm BUG end
194         ENDIF
195                                                            !
196         
197      ELSE                                             !* Leap-Frog : Asselin filter and swap
198         !                                ! =============!
199         IF( ln_linssh ) THEN             ! Fixed volume !
200            !                             ! =============!
201            DO jk = 1, jpkm1                             
202               DO jj = 1, jpj
203                  DO ji = 1, jpi   
204                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
205                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
206                     !
207                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
208                     vb(ji,jj,jk) = zvf
209                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
210                     vn(ji,jj,jk) = va(ji,jj,jk)
211                  END DO
212               END DO
213            END DO
214            !                             ! ================!
215         ELSE                             ! Variable volume !
216            !                             ! ================!
217            ! Before scale factor at t-points
218            ! (used as a now filtered scale factor until the swap)
219            ! ----------------------------------------------------
220            DO jk = 1, jpkm1
221               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
222            END DO
223            ! Add volume filter correction: compatibility with tracer advection scheme
224            ! => time filter + conservation correction (only at the first level)
225            zcoef = atfp * rdt * r1_rau0
226
227            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
228
229            IF ( ln_rnf ) THEN
230               IF( ln_rnf_depth ) THEN
231                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
232                     DO jj = 1, jpj
233                        DO ji = 1, jpi
234                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
235                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
236                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
237                           ENDIF
238                        ENDDO
239                     ENDDO
240                  ENDDO
241               ELSE
242                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
243               ENDIF
244            END IF
245
246            IF ( ln_isf ) THEN   ! if ice shelf melting
247               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
248                  DO jj = 1, jpj
249                     DO ji = 1, jpi
250                        IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN
251                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
252                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
253                        ELSEIF ( jk==misfkb(ji,jj) ) THEN
254                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
255                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
256                        ENDIF
257                     END DO
258                  END DO
259               END DO
260            END IF
261            !
262            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
263               ! Before filtered scale factor at (u/v)-points
264               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
265               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
266               DO jk = 1, jpkm1
267                  DO jj = 1, jpj
268                     DO ji = 1, jpi
269                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
270                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
271                        !
272                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
273                        vb(ji,jj,jk) = zvf
274                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
275                        vn(ji,jj,jk) = va(ji,jj,jk)
276                     END DO
277                  END DO
278               END DO
279               !
280            ELSE                          ! Asselin filter applied on thickness weighted velocity
281               !
282               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
283               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
284               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
285               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
286               DO jk = 1, jpkm1
287                  DO jj = 1, jpj
288                     DO ji = 1, jpi                 
289                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
290                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
291                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
292                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
293                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
294                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
295                        !
296                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
297                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
298                        !
299                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
300                        vb(ji,jj,jk) = zvf
301                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
302                        vn(ji,jj,jk) = va(ji,jj,jk)
303                     END DO
304                  END DO
305               END DO
306               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
307               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
308               !
309               DEALLOCATE( ze3u_f , ze3v_f )
310            ENDIF
311            !
312         ENDIF
313         !
314         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
315            ! Revert "before" velocities to time split estimate
316            ! Doing it here also means that asselin filter contribution is removed 
317            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
318            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
319            DO jk = 2, jpkm1
320               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
321               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
322            END DO
323            DO jk = 1, jpkm1
324               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
325               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
326            END DO
327         ENDIF
328         !
329      ENDIF ! neuler =/0
330      !
331      ! Set "now" and "before" barotropic velocities for next time step:
332      ! JC: Would be more clever to swap variables than to make a full vertical
333      ! integration
334      !
335      !
336      IF(.NOT.ln_linssh ) THEN
337         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
338         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
339         DO jk = 2, jpkm1
340            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
341            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
342         END DO
343         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
344         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
345      ENDIF
346      !
347      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
348      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
349      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
350      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
351      DO jk = 2, jpkm1
352         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
353         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
354         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
355         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
356      END DO
357      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
358      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
359      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
360      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
361      !
362      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
363         CALL iom_put(  "ubar", un_b(:,:) )
364         CALL iom_put(  "vbar", vn_b(:,:) )
365      ENDIF
366      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
367         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
368         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
369         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
370      ENDIF
371      !
372      IF ( iom_use("utau") ) THEN
373         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
374            ALLOCATE(zutau(jpi,jpj)) 
375            DO jj = 2, jpjm1
376               DO ji = 2, jpim1
377                  jk = miku(ji,jj) 
378                  zutau(ji,jj) = utau(ji,jj) & 
379                  &  + 0.5_wp * rau0 * (rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) * ua(ji,jj,jk) 
380               END DO
381            END DO
382            CALL lbc_lnk( 'dynnxt' , zutau, 'U', -1.)
383            CALL iom_put(  "utau", zutau(:,:) )
384            DEALLOCATE(zutau)
385         ELSE
386            CALL iom_put(  "utau", utau(:,:) )
387         ENDIF
388      ENDIF
389      !
390      IF ( iom_use("vtau") ) THEN
391         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
392            ALLOCATE(zvtau(jpi,jpj))
393            DO jj = 2, jpjm1
394               DO ji = 2, jpim1
395                  jk = mikv(ji,jj)
396                  zvtau(ji,jj) = vtau(ji,jj) &
397                  &  + 0.5_wp * rau0 * (rCdU_top(ji,jj+1)+rCdU_top(ji,jj)) * va(ji,jj,jk)
398               END DO
399            END DO
400            CALL lbc_lnk( 'dynnxt' , zvtau, 'V', -1.)
401            CALL iom_put(  "vtau", zvtau(:,:) )
402            DEALLOCATE(zvtau)
403         ELSE
404            CALL iom_put(  "vtau", vtau(:,:) )
405         ENDIF
406      ENDIF
407      !
408      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
409         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
410      !
411      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
412      IF( l_trddyn     )   DEALLOCATE( zua, zva )
413      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
414      !
415   END SUBROUTINE dyn_nxt
416
417   !!=========================================================================
418END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.