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 @ 12453

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

First implementation of the branch - compiling after merge
File size: 21.9 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   !!----------------------------------------------------------------------
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(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
103      !!----------------------------------------------------------------------
104      !
105      IF( ln_timing    )   CALL timing_start('dyn_nxt')
106      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
107      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
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(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
136               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + 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_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
148      !
149      !                                !* BDY open boundaries
150      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
151      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
152
153!!$   Do we need a call to bdy_vol here??
154      !
155      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
156         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
157         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
158         !
159         !                                  ! Kinetic energy and Conversion
160         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
161         !
162         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
163            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
164            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
165            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
166            CALL iom_put( "vtrd_tot", zva )
167         ENDIF
168         !
169         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
170         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
171         !                                  !  computation of the asselin filter trends)
172      ENDIF
173
174      ! Time filter and swap of dynamics arrays
175      ! ------------------------------------------
176      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
177         DO jk = 1, jpkm1
178            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
179            vn(:,:,jk) = va(:,:,jk)
180         END DO
181         ! limit velocities
182         IF (ln_ulimit) THEN
183            call dyn_limit_velocity (kt) 
184         ENDIF
185         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
186!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a 
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               !
192               e3t_n(:,:,jk) = e3t_a(:,:,jk)
193               e3u_n(:,:,jk) = e3u_a(:,:,jk)
194               e3v_n(:,:,jk) = e3v_a(:,:,jk)
195            END DO
196!!gm BUG end
197         ENDIF
198                                                            !
199         
200      ELSE                                             !* Leap-Frog : Asselin filter and swap
201         !                                ! =============!
202         IF( ln_linssh ) THEN             ! Fixed volume !
203            !                             ! =============!
204            DO jk = 1, jpkm1                             
205               DO jj = 1, jpj
206                  DO ji = 1, jpi   
207                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
208                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
209                     !
210                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
211                     vb(ji,jj,jk) = zvf
212                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
213                     vn(ji,jj,jk) = va(ji,jj,jk)
214                  END DO
215               END DO
216            END DO
217            ! limit velocities
218            IF (ln_ulimit) THEN
219               call dyn_limit_velocity (kt) 
220            ENDIF
221            !                             ! ================!
222         ELSE                             ! Variable volume !
223            !                             ! ================!
224            ! Before scale factor at t-points
225            ! (used as a now filtered scale factor until the swap)
226            ! ----------------------------------------------------
227            DO jk = 1, jpkm1
228               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
229            END DO
230            ! Add volume filter correction: compatibility with tracer advection scheme
231            ! => time filter + conservation correction (only at the first level)
232            zcoef = atfp * rdt * r1_rau0
233
234            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
235
236            IF ( ln_rnf ) THEN
237               IF( ln_rnf_depth ) THEN
238                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
239                     DO jj = 1, jpj
240                        DO ji = 1, jpi
241                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
242                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
243                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
244                           ENDIF
245                        ENDDO
246                     ENDDO
247                  ENDDO
248               ELSE
249                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
250               ENDIF
251            END IF
252
253            IF ( ln_isf ) THEN   ! if ice shelf melting
254               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
255                  DO jj = 1, jpj
256                     DO ji = 1, jpi
257                        IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN
258                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
259                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
260                        ELSEIF ( jk==misfkb(ji,jj) ) THEN
261                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
262                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
263                        ENDIF
264                     END DO
265                  END DO
266               END DO
267            END IF
268            !
269            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
270               ! Before filtered scale factor at (u/v)-points
271               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
272               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
273               DO jk = 1, jpkm1
274                  DO jj = 1, jpj
275                     DO ji = 1, jpi
276                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
277                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
278                        !
279                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
280                        vb(ji,jj,jk) = zvf
281                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
282                        vn(ji,jj,jk) = va(ji,jj,jk)
283                     END DO
284                  END DO
285               END DO
286               ! limit velocities
287               IF (ln_ulimit) THEN
288                  call dyn_limit_velocity (kt) 
289               ENDIF
290               !
291            ELSE                          ! Asselin filter applied on thickness weighted velocity
292               !
293               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
294               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
295               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
296               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
297               DO jk = 1, jpkm1
298                  DO jj = 1, jpj
299                     DO ji = 1, jpi                 
300                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
301                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
302                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
303                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
304                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
305                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
306                        !
307                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
308                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
309                        !
310                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
311                        vb(ji,jj,jk) = zvf
312                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
313                        vn(ji,jj,jk) = va(ji,jj,jk)
314                     END DO
315                  END DO
316               END DO
317               ! limit velocities
318               IF (ln_ulimit) THEN
319                  call dyn_limit_velocity (kt) 
320               ENDIF
321               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
322               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
323               !
324               DEALLOCATE( ze3u_f , ze3v_f )
325            ENDIF
326            !
327         ENDIF
328         !
329         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
330            ! Revert "before" velocities to time split estimate
331            ! Doing it here also means that asselin filter contribution is removed 
332            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
333            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
334            DO jk = 2, jpkm1
335               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
336               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
337            END DO
338            DO jk = 1, jpkm1
339               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
340               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
341            END DO
342         ENDIF
343         !
344      ENDIF ! neuler =/0
345      !
346      ! Set "now" and "before" barotropic velocities for next time step:
347      ! JC: Would be more clever to swap variables than to make a full vertical
348      ! integration
349      !
350      !
351      IF(.NOT.ln_linssh ) THEN
352         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
353         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
354         DO jk = 2, jpkm1
355            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
356            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
357         END DO
358         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
359         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
360      ENDIF
361      !
362      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
363      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
364      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
365      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
366      DO jk = 2, jpkm1
367         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
368         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
369         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
370         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
371      END DO
372      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
373      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
374      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
375      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
376      !
377      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
378         CALL iom_put(  "ubar", un_b(:,:) )
379         CALL iom_put(  "vbar", vn_b(:,:) )
380      ENDIF
381      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
382         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
383         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
384         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
385      ENDIF
386      !
387      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
388         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
389      !
390      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
391      IF( l_trddyn     )   DEALLOCATE( zua, zva )
392      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
393      !
394   END SUBROUTINE dyn_nxt
395
396   SUBROUTINE dyn_limit_velocity (kt) 
397      ! limits maxming vlaues of un and vn by volume courant number
398      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
399      !
400      INTEGER  ::   ji, jj, jk   ! dummy loop indices
401      REAL(wp) ::   zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm 
402 
403      ! limit fluxes
404      zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number
405      zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow
406 
407      DO jk = 1, jpkm1 
408         DO jj = 1, jpjm1 
409            DO ji = 1, jpim1 
410               ! U direction
411               zzu = un(ji,jj,jk) 
412               ze3e1 = e3u_n(ji  ,jj,jk) * e2u(ji  ,jj) 
413               ! ips is 1 if flow is positive othersize zero
414               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp ) 
415               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp ) 
416               ! idev is 1 if divergent flow otherwise zero
417               idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp ) 
418               idivm =  ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp ) 
419               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp 
420               zzcn = zzcn * zcn 
421               zplim =  zzcn * (e3t_n(ji  ,jj,jk) * e1t(ji  ,jj) * e2t(ji  ,jj)) / & 
422                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk) 
423               zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / & 
424                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk) 
425               ! limit currents
426               un(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism 
427               ! V  direction
428               zzu = vn(ji,jj,jk) 
429               ze3e1 = e3v_n(ji  ,jj,jk) * e1v(ji  ,jj) 
430               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp ) 
431               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp ) 
432               ! idev is 1 if divergent flow otherwise zero
433               idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp ) 
434               idivm =  ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp ) 
435               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp 
436               zzcn = zzcn * zcn 
437               zplim =  zzcn * (e3t_n(ji,jj  ,jk) * e1t(ji,jj  ) * e2t(ji,jj  )) / & 
438                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk) 
439               zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / & 
440                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk) 
441               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism 
442            ENDDO 
443         ENDDO 
444      ENDDO 
445 
446    END SUBROUTINE dyn_limit_velocity
447
448   !!=========================================================================
449END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.