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/r14075_India_uncoupled/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynnxt.F90 @ 15422

Last change on this file since 15422 was 15422, checked in by jcastill, 12 months ago

Changes tested so that they can merged with the CO9 Met Office branch - jpmax_harmo should be 34 with FES14 tides, but the last components are not used anyway

File size: 22.7 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 dynspg
36   USE domvvl         ! variable volume
37   USE bdy_oce , ONLY : ln_bdy
38   USE bdydta         ! ocean open boundary conditions
39   USE bdydyn         ! ocean open boundary conditions
40   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
41   USE trd_oce        ! trends: ocean variables
42   USE trddyn         ! trend manager: dynamics
43   USE trdken         ! trend manager: kinetic energy
44   !
45   USE in_out_manager ! I/O manager
46   USE iom            ! I/O manager library
47   USE lbclnk         ! lateral boundary condition (or mpp link)
48   USE lib_mpp        ! MPP library
49   USE prtctl         ! Print control
50   USE timing         ! Timing
51   USE zdfdrg ,  ONLY : ln_drgice_imp, rCdU_top
52#if defined key_agrif
53   USE agrif_oce_interp
54#endif
55
56   IMPLICIT NONE
57   PRIVATE
58
59   PUBLIC    dyn_nxt   ! routine called by step.F90
60
61   !! Substitution
62#  include "vectopt_loop_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
65   !! $Id$
66   !! Software governed by the CeCILL license (see ./LICENSE)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE dyn_nxt ( kt )
71      !!----------------------------------------------------------------------
72      !!                  ***  ROUTINE dyn_nxt  ***
73      !!                   
74      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
75      !!             condition on the after velocity, achieve the time stepping
76      !!             by applying the Asselin filter on now fields and swapping
77      !!             the fields.
78      !!
79      !! ** Method  : * Ensure after velocities transport matches time splitting
80      !!             estimate (ln_dynspg_ts=T)
81      !!
82      !!              * Apply lateral boundary conditions on after velocity
83      !!             at the local domain boundaries through lbc_lnk call,
84      !!             at the one-way open boundaries (ln_bdy=T),
85      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
86      !!
87      !!              * Apply the time filter applied and swap of the dynamics
88      !!             arrays to start the next time step:
89      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
90      !!                (un,vn) = (ua,va).
91      !!             Note that with flux form advection and non linear free surface,
92      !!             the time filter is applied on thickness weighted velocity.
93      !!             As a result, dyn_nxt MUST be called after tra_nxt.
94      !!
95      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
96      !!               un,vn   now horizontal velocity of next time-step
97      !!----------------------------------------------------------------------
98      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
99      !
100      INTEGER  ::   ji, jj, jk   ! dummy loop indices
101      INTEGER  ::   ikt          ! local integers
102      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
103      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
104      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
105      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zutau, zvtau
106      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
107      !!----------------------------------------------------------------------
108      !
109      IF( ln_timing    )   CALL timing_start('dyn_nxt')
110      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
111      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
112      !
113      IF( kt == nit000 ) THEN
114         IF(lwp) WRITE(numout,*)
115         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
116         IF(lwp) WRITE(numout,*) '~~~~~~~'
117      ENDIF
118
119      IF ( ln_dynspg_ts ) THEN
120         ! Ensure below that barotropic velocities match time splitting estimate
121         ! Compute actual transport and replace it with ts estimate at "after" time step
122         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
123         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
124         DO jk = 2, jpkm1
125            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
126            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
127         END DO
128         DO jk = 1, jpkm1
129            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
130            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
131         END DO
132         !
133         IF( .NOT.ln_bt_fw ) THEN
134            ! Remove advective velocity from "now velocities"
135            ! prior to asselin filtering     
136            ! In the forward case, this is done below after asselin filtering   
137            ! so that asselin contribution is removed at the same time
138            DO jk = 1, jpkm1
139               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
140               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
141            END DO 
142         ENDIF
143      ENDIF
144
145      ! Update after velocity on domain lateral boundaries
146      ! --------------------------------------------------     
147# if defined key_agrif
148      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
149# endif
150      !
151      CALL lbc_lnk_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
152      !
153      !                                !* BDY open boundaries
154      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
155      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
156
157!!$   Do we need a call to bdy_vol here??
158      !
159      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
160         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
161         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
162         !
163         !                                  ! Kinetic energy and Conversion
164         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
165         !
166         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
167            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
168            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
169            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
170            CALL iom_put( "vtrd_tot", zva )
171         ENDIF
172         !
173         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
174         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
175         !                                  !  computation of the asselin filter trends)
176      ENDIF
177
178      ! Time filter and swap of dynamics arrays
179      ! ------------------------------------------
180      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
181         DO jk = 1, jpkm1
182            ub(:,:,jk) = un(:,:,jk)                         ! ub <-- un
183            vb(:,:,jk) = vn(:,:,jk)
184            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
185            vn(:,:,jk) = va(:,:,jk)
186         END DO
187         ! limit velocities 
188         IF (ln_ulimit) THEN 
189            call dyn_limit_velocity (kt) 
190         ENDIF
191         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
192!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a 
193            DO jk = 1, jpkm1
194!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
195!               e3u_b(:,:,jk) = e3u_n(:,:,jk)
196!               e3v_b(:,:,jk) = e3v_n(:,:,jk)
197               !
198               e3t_n(:,:,jk) = e3t_a(:,:,jk)
199               e3u_n(:,:,jk) = e3u_a(:,:,jk)
200               e3v_n(:,:,jk) = e3v_a(:,:,jk)
201            END DO
202!!gm BUG end
203         ENDIF
204                                                            !
205         
206      ELSE                                             !* Leap-Frog : Asselin filter and swap
207         !                                ! =============!
208         IF( ln_linssh ) THEN             ! Fixed volume !
209            !                             ! =============!
210            DO jk = 1, jpkm1                             
211               DO jj = 1, jpj
212                  DO ji = 1, jpi   
213                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
214                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
215                     !
216                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
217                     vb(ji,jj,jk) = zvf
218                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
219                     vn(ji,jj,jk) = va(ji,jj,jk)
220                  END DO
221               END DO
222            END DO
223            ! limit velocities 
224            IF (ln_ulimit) THEN 
225               call dyn_limit_velocity (kt) 
226            ENDIF
227            !                             ! ================!
228         ELSE                             ! Variable volume !
229            !                             ! ================!
230            ! Before scale factor at t-points
231            ! (used as a now filtered scale factor until the swap)
232            ! ----------------------------------------------------
233            DO jk = 1, jpkm1
234               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
235            END DO
236            ! Add volume filter correction: compatibility with tracer advection scheme
237            ! => time filter + conservation correction (only at the first level)
238            zcoef = atfp * rdt * r1_rau0
239
240            DO jk = 1, jpkm1
241               e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,jk) & 
242                             &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
243            END DO
244
245            IF ( ln_rnf ) THEN
246               DO jk = 1, jpkm1
247                  e3t_b(:,:,jk) = e3t_b(:,:,jk) + zcoef * ( rnf_b(:,:) - rnf(:,:) ) * tmask(:,:,jk) & 
248                                &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
249               END DO
250            ENDIF
251
252            IF ( ln_isf ) THEN
253               DO jk = 1, jpkm1
254                  e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( fwfisf_b(:,:) - fwfisf(:,:) ) * tmask(:,:,jk) & 
255                                &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
256               END DO
257            ENDIF
258            !
259            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
260               ! Before filtered scale factor at (u/v)-points
261               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
262               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
263               DO jk = 1, jpkm1
264                  DO jj = 1, jpj
265                     DO ji = 1, jpi
266                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
267                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
268                        !
269                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
270                        vb(ji,jj,jk) = zvf
271                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
272                        vn(ji,jj,jk) = va(ji,jj,jk)
273                     END DO
274                  END DO
275               END DO
276               ! limit velocities 
277               IF (ln_ulimit) THEN 
278                  call dyn_limit_velocity (kt) 
279               ENDIF
280               !
281            ELSE                          ! Asselin filter applied on thickness weighted velocity
282               !
283               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
284               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
285               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
286               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
287               DO jk = 1, jpkm1
288                  DO jj = 1, jpj
289                     DO ji = 1, jpi                 
290                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
291                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
292                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
293                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
294                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
295                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
296                        !
297                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
298                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
299                        !
300                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
301                        vb(ji,jj,jk) = zvf
302                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
303                        vn(ji,jj,jk) = va(ji,jj,jk)
304                     END DO
305                  END DO
306               END DO
307               ! limit velocities 
308               IF (ln_ulimit) THEN 
309                  call dyn_limit_velocity (kt) 
310               ENDIF
311               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
312               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
313               !
314               DEALLOCATE( ze3u_f , ze3v_f )
315            ENDIF
316            !
317         ENDIF
318         !
319         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
320            ! Revert "before" velocities to time split estimate
321            ! Doing it here also means that asselin filter contribution is removed 
322            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
323            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
324            DO jk = 2, jpkm1
325               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
326               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
327            END DO
328            DO jk = 1, jpkm1
329               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
330               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
331            END DO
332         ENDIF
333         !
334      ENDIF ! neuler =/0
335      !
336      ! Set "now" and "before" barotropic velocities for next time step:
337      ! JC: Would be more clever to swap variables than to make a full vertical
338      ! integration
339      !
340      !
341      IF(.NOT.ln_linssh ) THEN
342         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
343         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
344         DO jk = 2, jpkm1
345            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
346            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
347         END DO
348         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
349         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
350      ENDIF
351      !
352      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
353      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
354      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
355      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
356      DO jk = 2, jpkm1
357         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
358         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
359         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
360         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
361      END DO
362      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
363      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
364      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
365      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
366      !
367      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
368         CALL iom_put(  "ubar", un_b(:,:) )
369         CALL iom_put(  "vbar", vn_b(:,:) )
370      ENDIF
371      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
372         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
373         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
374         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
375      ENDIF
376      !
377      IF ( iom_use("utau") ) THEN
378         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
379            ALLOCATE(zutau(jpi,jpj)) 
380            DO jj = 2, jpjm1
381               DO ji = 2, jpim1
382                  jk = miku(ji,jj) 
383                  zutau(ji,jj) = utau(ji,jj) & 
384                  &  + 0.5_wp * rau0 * (rCdU_top(ji+1,jj)+rCdU_top(ji,jj)) * ua(ji,jj,jk) 
385               END DO
386            END DO
387            CALL lbc_lnk( 'dynnxt' , zutau, 'U', -1.)
388            CALL iom_put(  "utau", zutau(:,:) )
389            DEALLOCATE(zutau)
390         ELSE
391            CALL iom_put(  "utau", utau(:,:) )
392         ENDIF
393      ENDIF
394      !
395      IF ( iom_use("vtau") ) THEN
396         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
397            ALLOCATE(zvtau(jpi,jpj))
398            DO jj = 2, jpjm1
399               DO ji = 2, jpim1
400                  jk = mikv(ji,jj)
401                  zvtau(ji,jj) = vtau(ji,jj) &
402                  &  + 0.5_wp * rau0 * (rCdU_top(ji,jj+1)+rCdU_top(ji,jj)) * va(ji,jj,jk)
403               END DO
404            END DO
405            CALL lbc_lnk( 'dynnxt' , zvtau, 'V', -1.)
406            CALL iom_put(  "vtau", zvtau(:,:) )
407            DEALLOCATE(zvtau)
408         ELSE
409            CALL iom_put(  "vtau", vtau(:,:) )
410         ENDIF
411      ENDIF
412      !
413      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
414         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
415      !
416      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
417      IF( l_trddyn     )   DEALLOCATE( zua, zva )
418      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
419      !
420   END SUBROUTINE dyn_nxt
421
422   SUBROUTINE dyn_limit_velocity (kt) 
423      ! limits maximum values of un and vn by volume courant number 
424      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
425     
426      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
427      REAL(wp) ::   zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm 
428   
429      ! limit fluxes 
430      zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number 
431      zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow 
432   
433      DO jk = 1, jpkm1 
434         DO jj = 2, jpjm1 
435            DO ji = fs_2, fs_jpim1   ! vect. opt.
436               ! U direction 
437               zzu = un(ji,jj,jk) 
438               ze3e1 = e3u_n(ji  ,jj,jk) * e2u(ji  ,jj) 
439               ! ips is 1 if flow is positive othersize zero 
440               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp ) 
441               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp ) 
442               ! idev is 1 if divergent flow otherwise zero 
443               idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp ) 
444               idivm =  ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp ) 
445               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp 
446               zzcn = zzcn * zcn 
447               zplim =  zzcn * (e3t_n(ji  ,jj,jk) * e1t(ji  ,jj) * e2t(ji  ,jj)) / & 
448                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk) 
449               zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / & 
450                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk) 
451               ! limit currents 
452               un(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism 
453   
454               ! V  direction 
455               zzu = vn(ji,jj,jk) 
456               ze3e1 = e3v_n(ji  ,jj,jk) * e1v(ji  ,jj) 
457               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp ) 
458               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp ) 
459               ! idev is 1 if divergent flow otherwise zero 
460               idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp ) 
461               idivm =  ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp ) 
462               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp 
463               zzcn = zzcn * zcn 
464               zplim =  zzcn * (e3t_n(ji,jj  ,jk) * e1t(ji,jj  ) * e2t(ji,jj  )) / & 
465                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk) 
466               zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / & 
467                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk) 
468               ! limit currents 
469               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism 
470            ENDDO 
471         ENDDO 
472      ENDDO 
473       CALL lbc_lnk_multi( 'dynnxt', un(:,:,:), 'U',  -1., vn(:,:,:), 'V',  -1. ) 
474   
475    END SUBROUTINE dyn_limit_velocity
476
477   !!=========================================================================
478END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.