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

source: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 11080

Last change on this file since 11080 was 11080, checked in by jcastill, 5 years ago

Latest changes from Jason
File size: 20.8 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 dynspg
34   USE domvvl         ! variable volume
35   USE bdy_oce   , ONLY: ln_bdy
36   USE bdydta         ! ocean open boundary conditions
37   USE bdydyn         ! ocean open boundary conditions
38   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
39   USE trd_oce        ! trends: ocean variables
40   USE trddyn         ! trend manager: dynamics
41   USE trdken         ! trend manager: kinetic energy
42   !
43   USE in_out_manager ! I/O manager
44   USE iom            ! I/O manager library
45   USE lbclnk         ! lateral boundary condition (or mpp link)
46   USE lib_mpp        ! MPP library
47   USE wrk_nemo       ! Memory Allocation
48   USE prtctl         ! Print control
49   USE timing         ! Timing
50#if defined key_agrif
51   USE agrif_opa_interp
52#endif
53
54   IMPLICIT NONE
55   PRIVATE
56
57   PUBLIC    dyn_nxt   ! routine called by step.F90
58
59   !!----------------------------------------------------------------------
60   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE dyn_nxt ( kt )
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      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
86      !!                (un,vn) = (ua,va).
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 :   ub,vb   filtered before horizontal velocity of next time-step
92      !!               un,vn   now horizontal velocity of next time-step
93      !!----------------------------------------------------------------------
94      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
95      !
96      INTEGER  ::   ji, jj, jk   ! dummy loop indices
97      INTEGER  ::   ikt          ! local integers
98      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
99      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
100      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
101      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
102      !!----------------------------------------------------------------------
103      !
104      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
105      !
106      IF( ln_dynspg_ts       )   CALL wrk_alloc( jpi,jpj,       zue, zve)
107      IF( l_trddyn           )   CALL wrk_alloc( jpi,jpj,jpk,   zua, zva)
108      !
109      IF( kt == nit000 ) THEN
110         IF(lwp) WRITE(numout,*)
111         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
112         IF(lwp) WRITE(numout,*) '~~~~~~~'
113      ENDIF
114
115      IF ( ln_dynspg_ts ) THEN
116         ! Ensure below that barotropic velocities match time splitting estimate
117         ! Compute actual transport and replace it with ts estimate at "after" time step
118         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
119         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
120         DO jk = 2, jpkm1
121            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
122            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
123         END DO
124         DO jk = 1, jpkm1
125            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
126            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
127         END DO
128         !
129         IF( .NOT.ln_bt_fw ) THEN
130            ! Remove advective velocity from "now velocities"
131            ! prior to asselin filtering     
132            ! In the forward case, this is done below after asselin filtering   
133            ! so that asselin contribution is removed at the same time
134            DO jk = 1, jpkm1
135               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
136               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
137            END DO 
138         ENDIF
139      ENDIF
140
141      ! Update after velocity on domain lateral boundaries
142      ! --------------------------------------------------     
143# if defined key_agrif
144      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
145# endif
146      !
147      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
148      CALL lbc_lnk( va, 'V', -1. ) 
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         ! limit velocities
183         IF (ln_ulimit) THEN
184            call dyn_limit_velocity (kt)
185         ENDIF
186         IF(.NOT.ln_linssh ) THEN
187            DO jk = 1, jpkm1
188               e3t_b(:,:,jk) = e3t_n(:,:,jk)
189               e3u_b(:,:,jk) = e3u_n(:,:,jk)
190               e3v_b(:,:,jk) = e3v_n(:,:,jk)
191            END DO
192         ENDIF
193      ELSE                                             !* Leap-Frog : Asselin filter and swap
194         !                                ! =============!
195         IF( ln_linssh ) THEN             ! Fixed volume !
196            !                             ! =============!
197            DO jk = 1, jpkm1                             
198               DO jj = 1, jpj
199                  DO ji = 1, jpi   
200                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
201                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
202                     !
203                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
204                     vb(ji,jj,jk) = zvf
205                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
206                     vn(ji,jj,jk) = va(ji,jj,jk)
207                  END DO
208               END DO
209            END DO
210            ! limit velocities
211            IF (ln_ulimit) THEN
212               call dyn_limit_velocity (kt)
213            ENDIF
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            IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN    ! No asselin filtering on thicknesses if forward time splitting
221               e3t_b(:,:,1:jpkm1) = e3t_n(:,:,1:jpkm1)
222            ELSE
223               DO jk = 1, jpkm1
224                  e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
225               END DO
226               ! Add volume filter correction: compatibility with tracer advection scheme
227               ! => time filter + conservation correction (only at the first level)
228               zcoef = atfp * rdt * r1_rau0
229               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting
230                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) &
231                                 &                      - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
232               ELSE                     ! if ice shelf melting
233                  DO jj = 1, jpj
234                     DO ji = 1, jpi
235                        ikt = mikt(ji,jj)
236                        e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  &
237                           &                                           - rnf_b   (ji,jj) + rnf   (ji,jj)  &
238                           &                                           + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,ikt)
239                     END DO
240                  END DO
241               END IF
242            ENDIF
243            !
244            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
245               ! Before filtered scale factor at (u/v)-points
246               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
247               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
248               DO jk = 1, jpkm1
249                  DO jj = 1, jpj
250                     DO ji = 1, jpi
251                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
252                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
253                        !
254                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
255                        vb(ji,jj,jk) = zvf
256                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
257                        vn(ji,jj,jk) = va(ji,jj,jk)
258                     END DO
259                  END DO
260               END DO
261               ! limit velocities
262               IF (ln_ulimit) THEN
263                  call dyn_limit_velocity (kt)
264               ENDIF
265               !
266            ELSE                          ! Asselin filter applied on thickness weighted velocity
267               !
268               CALL wrk_alloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
269               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
270               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
271               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
272               DO jk = 1, jpkm1
273                  DO jj = 1, jpj
274                     DO ji = 1, jpi                 
275                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
276                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
277                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
278                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
279                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
280                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
281                        !
282                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
283                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
284                        !
285                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
286                        vb(ji,jj,jk) = zvf
287                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
288                        vn(ji,jj,jk) = va(ji,jj,jk)
289                     END DO
290                  END DO
291               END DO
292               ! limit velocities
293               IF (ln_ulimit) THEN
294                  call dyn_limit_velocity (kt)
295               ENDIF
296               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
297               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
298               !
299               CALL wrk_dealloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
300            ENDIF
301            !
302         ENDIF
303         !
304         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
305            ! Revert "before" velocities to time split estimate
306            ! Doing it here also means that asselin filter contribution is removed 
307            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
308            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
309            DO jk = 2, jpkm1
310               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
311               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
312            END DO
313            DO jk = 1, jpkm1
314               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
315               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
316            END DO
317         ENDIF
318         !
319      ENDIF ! neuler =/0
320      !
321      ! Set "now" and "before" barotropic velocities for next time step:
322      ! JC: Would be more clever to swap variables than to make a full vertical
323      ! integration
324      !
325      !
326      IF(.NOT.ln_linssh ) THEN
327         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
328         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
329         DO jk = 2, jpkm1
330            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
331            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
332         END DO
333         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
334         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
335      ENDIF
336      !
337      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
338      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
339      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
340      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
341      DO jk = 2, jpkm1
342         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
343         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
344         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
345         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
346      END DO
347      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
348      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
349      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
350      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
351      !
352      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
353         CALL iom_put(  "ubar", un_b(:,:) )
354         CALL iom_put(  "vbar", vn_b(:,:) )
355      ENDIF
356      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
357         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
358         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
359         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
360      ENDIF
361      !
362      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
363         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
364      !
365      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj,       zue, zve )
366      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk,   zua, zva )
367      !
368      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
369      !
370   END SUBROUTINE dyn_nxt
371
372   SUBROUTINE dyn_limit_velocity (kt)
373      ! limits maxming vlaues of un and vn by volume courant number
374      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
375      !
376      INTEGER  ::   ji, jj, jk   ! dummy loop indices
377      REAL(wp) :: zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm
378
379      ! limit fluxes
380      zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number
381      zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow
382
383      DO jk = 1, jpkm1
384         DO jj = 1, jpjm1
385            DO ji = 1, jpim1
386               ! U direction
387               zzu = un(ji,jj,jk)
388               ze3e1 = e3u_n(ji  ,jj,jk) * e2u(ji  ,jj)
389               ! ips is 1 if flow is positive othersize zero
390               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )
391               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )
392               ! idev is 1 if divergent flow otherwise zero
393               idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp )
394               idivm =  ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp )
395               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp
396               zzcn = zzcn * zcn
397               zplim =  zzcn * (e3t_n(ji  ,jj,jk) * e1t(ji  ,jj) * e2t(ji  ,jj)) / &
398                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)
399               zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / &
400                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)
401               ! limit currents
402               un(ji,jj,jk) = min ( zzu,zplim) * isp + max (zzu,zmlim) *ism
403               ! V  direction
404               zzu = vn(ji,jj,jk)
405               ze3e1 = e3v_n(ji  ,jj,jk) * e1v(ji  ,jj)
406               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )
407               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )
408               ! idev is 1 if divergent flow otherwise zero
409               idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp )
410               idivm =  ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp )
411               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp
412               zzcn = zzcn * zcn
413               zplim =  zzcn * (e3t_n(ji,jj  ,jk) * e1t(ji,jj  ) * e2t(ji,jj  )) / &
414                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)
415               zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / &
416                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)
417               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max (zzu,zmlim) *ism
418            ENDDO
419         ENDDO
420      ENDDO
421
422    END SUBROUTINE dyn_limit_velocity
423
424   !!=========================================================================
425END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.