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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 4455

Last change on this file since 4455 was 4455, checked in by trackstand2, 10 years ago

Add mbkmax to dyn_nxt
  • Property svn:keywords set to Id
File size: 21.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   !!-------------------------------------------------------------------------
20 
21   !!-------------------------------------------------------------------------
22   !!   dyn_nxt      : obtain the next (after) horizontal velocity
23   !!-------------------------------------------------------------------------
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
26   USE sbc_oce         ! Surface boundary condition: ocean fields
27   USE phycst          ! physical constants
28   USE dynspg_oce      ! type of surface pressure gradient
29   USE dynadv          ! dynamics: vector invariant versus flux form
30   USE domvvl          ! variable volume
31   USE obc_oce         ! ocean open boundary conditions
32   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
33   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
34   USE obcvol          ! ocean open boundary condition (obc_vol routines)
35   USE bdy_oce         ! unstructured open boundary conditions
36   USE bdydta          ! unstructured open boundary conditions
37   USE bdydyn          ! unstructured open boundary conditions
38   USE in_out_manager  ! I/O manager
39   USE lbclnk          ! lateral boundary condition (or mpp link)
40   USE lib_mpp         ! MPP library
41   USE prtctl          ! Print control
42#if defined key_agrif
43   USE agrif_opa_interp
44#endif
45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC    dyn_nxt   ! routine called by step.F90
50
51   !! * Control permutation of array indices
52#  include "oce_ftrans.h90"
53#  include "dom_oce_ftrans.h90"
54#  include "sbc_oce_ftrans.h90"
55#  include "domvvl_ftrans.h90"
56#  include "obc_oce_ftrans.h90"
57
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60   !!----------------------------------------------------------------------
61   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
62   !! $Id$
63   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE dyn_nxt ( kt )
68      !!----------------------------------------------------------------------
69      !!                  ***  ROUTINE dyn_nxt  ***
70      !!                   
71      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
72      !!             condition on the after velocity, achieved the time stepping
73      !!             by applying the Asselin filter on now fields and swapping
74      !!             the fields.
75      !!
76      !! ** Method  : * After velocity is compute using a leap-frog scheme:
77      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
78      !!             Note that with flux form advection and variable volume layer
79      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
80      !!             velocity.
81      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
82      !!             the time stepping has already been done in dynspg module
83      !!
84      !!              * Apply lateral boundary conditions on after velocity
85      !!             at the local domain boundaries through lbc_lnk call,
86      !!             at the radiative open boundaries (lk_obc=T),
87      !!             at the relaxed   open boundaries (lk_bdy=T), and
88      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
89      !!
90      !!              * Apply the time filter applied and swap of the dynamics
91      !!             arrays to start the next time step:
92      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
93      !!                (un,vn) = (ua,va).
94      !!             Note that with flux form advection and variable volume layer
95      !!             (lk_vvl=T), the time filter is applied on thickness weighted
96      !!             velocity.
97      !!
98      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
99      !!               un,vn   now horizontal velocity of next time-step
100      !!----------------------------------------------------------------------
101      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
102      USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
103      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
104      USE arpdebugging, ONLY: dump_array
105      !! DCSE_NEMO: need additional directives for renamed module variables
106!FTRANS ze3u_f :I :I :z
107!FTRANS ze3v_f :I :I :z
108      !
109      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
110      !
111      INTEGER  ::   ji, jj, jk   ! dummy loop indices
112#if ! defined key_dynspg_flt
113      REAL(wp) ::   z2dt         ! temporary scalar
114#endif
115      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
116      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
117      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
118      !!----------------------------------------------------------------------
119
120      IF( wrk_in_use(2, 1,2,3) ) THEN
121         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
122      ENDIF
123
124      IF( kt == nit000 ) THEN
125         IF(lwp) WRITE(numout,*)
126         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
127         IF(lwp) WRITE(numout,*) '~~~~~~~'
128      ENDIF
129
130!      CALL dump_array(kt, 'ua_nxt_start',ua(:,:,1),withHalos=.TRUE.)
131
132#if defined key_dynspg_flt
133      !
134      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
135      ! -------------
136
137      ! Update after velocity on domain lateral boundaries      (only local domain required)
138      ! --------------------------------------------------
139      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
140      CALL lbc_lnk( va, 'V', -1. ) 
141      !
142#else
143      ! Next velocity :   Leap-frog time stepping
144      ! -------------
145      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
146      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
147      !
148      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
149         DO jk = 1, jpkm1
150            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
151            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
152         END DO
153      ELSE                                            ! applied on thickness weighted velocity
154
155#if defined key_z_first
156         DO jj = 1, jpj, 1
157            DO ji = 1, jpi, 1
158               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
159                  ua(ji,jj,jk) = ( ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      &
160               &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   &
161               &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
162                  va(ji,jj,jk) = ( vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      &
163               &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   &
164               &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
165               END DO
166            END DO
167         END DO
168#else
169         DO jk = 1, jpkm1
170            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
171               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
172               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
173            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
174               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
175               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
176         END DO
177#endif
178      ENDIF
179
180
181      ! Update after velocity on domain lateral boundaries
182      ! --------------------------------------------------     
183      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
184      CALL lbc_lnk( va, 'V', -1. ) 
185      !
186#if defined ARPDBGSUM
187      WRITE (*,*) 'ARPDBG: dyn_nxt: after lbc update sum(ua)=',SUM(ua),' at step=',kt
188#endif
189
190# if defined key_obc
191      !                                !* OBC open boundaries
192      CALL obc_dyn( kt )
193      !
194      IF( .NOT. lk_dynspg_flt ) THEN
195         ! Flather boundary condition : - Update sea surface height on each open boundary
196         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
197         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
198         !                              - Correct the barotropic velocities
199         CALL obc_dyn_bt( kt )
200         !
201!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
202         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
203         !
204         IF( ln_vol_cst )   CALL obc_vol( kt )
205         !
206         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
207      ENDIF
208      !
209# elif defined key_bdy 
210      !                                !* BDY open boundaries
211      IF( .NOT. lk_dynspg_flt ) THEN
212         CALL bdy_dyn_frs( kt )
213#  if ! defined key_vvl
214         ua_e(:,:) = 0.e0
215         va_e(:,:) = 0.e0
216         ! Set these variables for use in bdy_dyn_fla
217         hur_e(:,:) = hur(:,:)
218         hvr_e(:,:) = hvr(:,:)
219         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
220            ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
221            va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
222         END DO
223         ua_e(:,:) = ua_e(:,:) * hur(:,:)
224         va_e(:,:) = va_e(:,:) * hvr(:,:)
225         DO jk = 1 , jpkm1
226            ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
227            va(:,:,jk) = va(:,:,jk) - va_e(:,:)
228         END DO
229         CALL bdy_dta_fla( kt+1, 0,2*nn_baro)
230         CALL bdy_dyn_fla( sshn_b )
231         CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
232         CALL lbc_lnk( va_e, 'V', -1. )   !
233         DO jk = 1 , jpkm1
234            ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
235            va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
236         END DO
237#  endif
238      ENDIF
239# endif
240      !
241# if defined key_agrif
242      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
243# endif
244#endif
245
246#if defined ARPDBGSUM
247      WRITE(*,*) 'ARPDBG: dyn_nxt b4 time filter SUM(ua)=',SUM(ua),'at step=',kt
248#endif
249
250      ! Time filter and swap of dynamics arrays
251      ! ------------------------------------------
252      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
253#if defined key_z_first
254         DO jj = 1, jpj
255            DO ji = 1, jpi
256               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
257                  un(ji,jj,jk) = ua(ji,jj,jk)                ! un <-- ua
258                  vn(ji,jj,jk) = va(ji,jj,jk)
259               END DO
260            END DO
261         END DO
262#else
263         DO jk = 1, jpkm1
264            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
265            vn(:,:,jk) = va(:,:,jk)
266         END DO
267#endif
268      ELSE                                             !* Leap-Frog : Asselin filter and swap
269         !                                ! =============!
270         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
271            !                             ! =============!
272#if defined key_z_first
273            DO jj = 1, jpj
274               DO ji = 1, jpi   
275                  DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1                             
276#else
277            DO jk = 1, jpkm1                             
278               DO jj = 1, jpj
279                  DO ji = 1, jpi   
280#endif
281                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
282                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
283                     !
284                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
285                     vb(ji,jj,jk) = zvf
286                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
287                     vn(ji,jj,jk) = va(ji,jj,jk)
288                  END DO
289               END DO
290            END DO
291            !                             ! ================!
292         ELSE                             ! Variable volume !
293            !                             ! ================!
294            ! Before scale factor at t-points
295            ! -------------------------------
296#if defined key_z_first
297            DO jj = 1, jpj
298               DO ji = 1, jpi   
299                  DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
300                     fse3t_b(ji,jj,jk) = fse3t_n(ji,jj,jk)                             &
301                  &              + atfp * (  fse3t_b(ji,jj,jk) + fse3t_a(ji,jj,jk)     &
302                  &                         - 2.e0 * fse3t_n(ji,jj,jk)            )
303                  END DO
304               END DO
305            ENDDO
306#else
307            DO jk = 1, jpkm1
308               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
309                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
310                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
311            ENDDO
312#endif
313            ! Add volume filter correction only at the first level of t-point scale factors
314            zec = atfp * rdt / rau0
315#if defined key_z_first
316            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask_1(:,:)
317#else
318            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
319#endif
320            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
321            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
322            zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
323            zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
324            !
325            IF( ln_dynadv_vec ) THEN
326               ! Before scale factor at (u/v)-points
327               ! -----------------------------------
328               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
329#if defined key_z_first
330               DO jj = 1, jpjm1
331                  DO ji = 1, jpim1
332                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
333#else
334               DO jk = 1, jpkm1
335                  DO jj = 1, jpjm1
336                     DO ji = 1, jpim1
337#endif
338                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
339                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
340                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
341                        fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
342                        fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
343                     END DO
344                  END DO
345               END DO
346               ! lateral boundary conditions
347               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
348               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
349               ! Add initial scale factor to scale factor anomaly
350#if defined key_z_first
351               DO jj = 1, jpjm1
352                  DO ji = 1, jpim1
353                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
354                        fse3u_b(ji,jj,jk) = fse3u_b(ji,jj,jk) + fse3u_0(ji,jj,jk)
355                        fse3v_b(ji,jj,jk) = fse3v_b(ji,jj,jk) + fse3v_0(ji,jj,jk)
356                     END DO
357                  END DO
358               END DO
359#else
360               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
361               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
362#endif
363               ! Leap-Frog - Asselin filter and swap: applied on velocity
364               ! -----------------------------------
365#if defined key_z_first
366               DO jj = 1, jpj
367                  DO ji = 1, jpi
368                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
369#else
370               DO jk = 1, jpkm1
371                  DO jj = 1, jpj
372                     DO ji = 1, jpi
373#endif
374                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
375                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
376                        !
377                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
378                        vb(ji,jj,jk) = zvf
379                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
380                        vn(ji,jj,jk) = va(ji,jj,jk)
381                     END DO
382                  END DO
383               END DO
384               !
385            ELSE
386               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
387               !-----------------------------------------------
388               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
389#if defined key_z_first
390               DO jj = 1, jpjm1
391                  DO ji = 1, jpim1
392                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
393#else
394               DO jk = 1, jpkm1
395                  DO jj = 1, jpjm1
396                     DO ji = 1, jpim1
397#endif
398                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
399                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
400                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
401                        ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) )
402                        ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) )
403                     END DO
404                  END DO
405               END DO
406               ! lateral boundary conditions
407               CALL lbc_lnk( ze3u_f, 'U', 1. )
408               CALL lbc_lnk( ze3v_f, 'V', 1. )
409               ! Add initial scale factor to scale factor anomaly
410#if defined key_z_first
411               DO jj = 1, jpjm1
412                  DO ji = 1, jpim1
413                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
414                        ze3u_f(ji,jj,jk) = ze3u_f(ji,jj,jk) + fse3u_0(ji,jj,jk)
415                        ze3v_f(ji,jj,jk) = ze3v_f(ji,jj,jk) + fse3v_0(ji,jj,jk)
416                     END DO
417                  END DO
418               END DO
419#else
420               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
421               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
422#endif
423               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
424               ! -----------------------------------             ===========================
425#if defined key_z_first
426               DO jj = 1, jpj
427                  DO ji = 1, jpim1
428                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
429#else
430               DO jk = 1, jpkm1
431                  DO jj = 1, jpj
432                     DO ji = 1, jpim1
433#endif
434                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
435                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
436                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
437                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
438                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
439                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
440                        !
441                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
442                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
443                        !
444                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
445                        vb(ji,jj,jk) = zvf
446                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
447                        vn(ji,jj,jk) = va(ji,jj,jk)
448                     END DO
449                  END DO
450               END DO
451#if defined key_z_first
452               DO jj = 1, jpj
453                  DO ji = 1, jpim1
454                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
455                        fse3u_b(ji,jj,jk) = ze3u_f(ji,jj,jk)                   ! e3u_b <-- filtered scale factor
456                        fse3v_b(ji,jj,jk) = ze3v_f(ji,jj,jk)
457                     END DO
458                  END DO
459               END DO
460#else
461               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
462               fse3v_b(:,:,:) = ze3v_f(:,:,:)
463#endif
464               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
465               CALL lbc_lnk( vb, 'V', -1. )
466            ENDIF
467            !
468         ENDIF
469         !
470      ENDIF
471
472#if defined ARPDBGSUM
473      WRITE(*,*) 'ARPDBG: dyn_nxt end SUM(un)=',SUM(un),'at step=',kt
474#endif
475
476      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
477         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
478      !
479      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
480      !
481   END SUBROUTINE dyn_nxt
482
483   !!=========================================================================
484END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.