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

source: NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DYN/dynnxt.F90 @ 12554

Last change on this file since 12554 was 12554, checked in by jcastill, 4 years ago

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