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

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

Add debug SUM output to dynnxt and Dynamics section of step

  • Property svn:keywords set to Id
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   !!-------------------------------------------------------------------------
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#if defined ARPDBGSUM
132      WRITE (*,*) narea,': ARPDBG: dyn_nxt: start, sum(ua)=',SUM(ua),' at step=',kt
133#endif
134
135#if defined key_dynspg_flt
136      !
137      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
138      ! -------------
139
140      ! Update after velocity on domain lateral boundaries      (only local domain required)
141      ! --------------------------------------------------
142      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
143      CALL lbc_lnk( va, 'V', -1. ) 
144      !
145#else
146      ! Next velocity :   Leap-frog time stepping
147      ! -------------
148      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
149      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
150      !
151      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
152         DO jk = 1, jpkm1
153            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
154            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
155         END DO
156      ELSE                                            ! applied on thickness weighted velocity
157
158#if defined key_z_first
159         DO jj = 1, jpj, 1
160            DO ji = 1, jpi, 1
161               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
162                  ua(ji,jj,jk) = ( ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      &
163               &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   &
164               &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
165                  va(ji,jj,jk) = ( vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      &
166               &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   &
167               &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
168               END DO
169            END DO
170         END DO
171#else
172         DO jk = 1, jpkm1
173            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
174               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
175               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
176            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
177               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
178               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
179         END DO
180#endif
181      ENDIF
182
183
184      ! Update after velocity on domain lateral boundaries
185      ! --------------------------------------------------     
186      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
187      CALL lbc_lnk( va, 'V', -1. ) 
188      !
189#if defined ARPDBGSUM
190      WRITE (*,*) narea,': ARPDBG: dyn_nxt: after lbc update sum(ua)=',SUM(ua),' at step=',kt
191#endif
192
193# if defined key_obc
194      !                                !* OBC open boundaries
195      CALL obc_dyn( kt )
196      !
197      IF( .NOT. lk_dynspg_flt ) THEN
198         ! Flather boundary condition : - Update sea surface height on each open boundary
199         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
200         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
201         !                              - Correct the barotropic velocities
202         CALL obc_dyn_bt( kt )
203         !
204!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
205         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
206         !
207         IF( ln_vol_cst )   CALL obc_vol( kt )
208         !
209         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
210      ENDIF
211      !
212# elif defined key_bdy 
213      !                                !* BDY open boundaries
214      IF( .NOT. lk_dynspg_flt ) THEN
215         CALL bdy_dyn_frs( kt )
216#  if ! defined key_vvl
217         ua_e(:,:) = 0.e0
218         va_e(:,:) = 0.e0
219         ! Set these variables for use in bdy_dyn_fla
220         hur_e(:,:) = hur(:,:)
221         hvr_e(:,:) = hvr(:,:)
222         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
223            ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
224            va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
225         END DO
226         ua_e(:,:) = ua_e(:,:) * hur(:,:)
227         va_e(:,:) = va_e(:,:) * hvr(:,:)
228         DO jk = 1 , jpkm1
229            ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
230            va(:,:,jk) = va(:,:,jk) - va_e(:,:)
231         END DO
232         CALL bdy_dta_fla( kt+1, 0,2*nn_baro)
233         CALL bdy_dyn_fla( sshn_b )
234         CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
235         CALL lbc_lnk( va_e, 'V', -1. )   !
236         DO jk = 1 , jpkm1
237            ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
238            va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
239         END DO
240#  endif
241      ENDIF
242# endif
243      !
244# if defined key_agrif
245      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
246# endif
247#endif
248
249#if defined ARPDBGSUM
250      WRITE(*,*) narea,': ARPDBG: dyn_nxt b4 time filter SUM(ua)=',SUM(ua),'at step=',kt
251#endif
252
253      ! Time filter and swap of dynamics arrays
254      ! ------------------------------------------
255      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
256#if defined key_z_first
257         DO jj = 1, jpj
258            DO ji = 1, jpi
259               DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
260                  un(ji,jj,jk) = ua(ji,jj,jk)                ! un <-- ua
261                  vn(ji,jj,jk) = va(ji,jj,jk)
262               END DO
263            END DO
264         END DO
265#else
266         DO jk = 1, jpkm1
267            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
268            vn(:,:,jk) = va(:,:,jk)
269         END DO
270#endif
271      ELSE                                             !* Leap-Frog : Asselin filter and swap
272         !                                ! =============!
273         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
274            !                             ! =============!
275#if defined key_z_first
276            DO jj = 1, jpj
277               DO ji = 1, jpi   
278                  DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1                             
279#else
280            DO jk = 1, jpkm1                             
281               DO jj = 1, jpj
282                  DO ji = 1, jpi   
283#endif
284                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
285                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
286                     !
287                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
288                     vb(ji,jj,jk) = zvf
289                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
290                     vn(ji,jj,jk) = va(ji,jj,jk)
291                  END DO
292               END DO
293            END DO
294            !                             ! ================!
295         ELSE                             ! Variable volume !
296            !                             ! ================!
297            ! Before scale factor at t-points
298            ! -------------------------------
299#if defined key_z_first
300            DO jj = 1, jpj
301               DO ji = 1, jpi   
302                  DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
303                     fse3t_b(ji,jj,jk) = fse3t_n(ji,jj,jk)                             &
304                  &              + atfp * (  fse3t_b(ji,jj,jk) + fse3t_a(ji,jj,jk)     &
305                  &                         - 2.e0 * fse3t_n(ji,jj,jk)            )
306                  END DO
307               END DO
308            ENDDO
309#else
310            DO jk = 1, jpkm1
311               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
312                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
313                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
314            ENDDO
315#endif
316            ! Add volume filter correction only at the first level of t-point scale factors
317            zec = atfp * rdt / rau0
318#if defined key_z_first
319            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask_1(:,:)
320#else
321            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
322#endif
323            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
324            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
325            zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
326            zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
327            !
328            IF( ln_dynadv_vec ) THEN
329               ! Before scale factor at (u/v)-points
330               ! -----------------------------------
331               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
332#if defined key_z_first
333               DO jj = 1, jpjm1
334                  DO ji = 1, jpim1
335                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
336#else
337               DO jk = 1, jpkm1
338                  DO jj = 1, jpjm1
339                     DO ji = 1, jpim1
340#endif
341                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
342                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
343                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
344                        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) )
345                        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) )
346                     END DO
347                  END DO
348               END DO
349               ! lateral boundary conditions
350               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
351               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
352               ! Add initial scale factor to scale factor anomaly
353#if defined key_z_first
354               DO jj = 1, jpjm1
355                  DO ji = 1, jpim1
356                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
357                        fse3u_b(ji,jj,jk) = fse3u_b(ji,jj,jk) + fse3u_0(ji,jj,jk)
358                        fse3v_b(ji,jj,jk) = fse3v_b(ji,jj,jk) + fse3v_0(ji,jj,jk)
359                     END DO
360                  END DO
361               END DO
362#else
363               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
364               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
365#endif
366               ! Leap-Frog - Asselin filter and swap: applied on velocity
367               ! -----------------------------------
368#if defined key_z_first
369               DO jj = 1, jpj
370                  DO ji = 1, jpi
371                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
372#else
373               DO jk = 1, jpkm1
374                  DO jj = 1, jpj
375                     DO ji = 1, jpi
376#endif
377                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
378                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
379                        !
380                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
381                        vb(ji,jj,jk) = zvf
382                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
383                        vn(ji,jj,jk) = va(ji,jj,jk)
384                     END DO
385                  END DO
386               END DO
387               !
388            ELSE
389               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
390               !-----------------------------------------------
391               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
392#if defined key_z_first
393               DO jj = 1, jpjm1
394                  DO ji = 1, jpim1
395                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
396#else
397               DO jk = 1, jpkm1
398                  DO jj = 1, jpjm1
399                     DO ji = 1, jpim1
400#endif
401                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
402                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
403                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
404                        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) )
405                        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) )
406                     END DO
407                  END DO
408               END DO
409               ! lateral boundary conditions
410               CALL lbc_lnk( ze3u_f, 'U', 1. )
411               CALL lbc_lnk( ze3v_f, 'V', 1. )
412               ! Add initial scale factor to scale factor anomaly
413#if defined key_z_first
414               DO jj = 1, jpjm1
415                  DO ji = 1, jpim1
416                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
417                        ze3u_f(ji,jj,jk) = ze3u_f(ji,jj,jk) + fse3u_0(ji,jj,jk)
418                        ze3v_f(ji,jj,jk) = ze3v_f(ji,jj,jk) + fse3v_0(ji,jj,jk)
419                     END DO
420                  END DO
421               END DO
422#else
423               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
424               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
425#endif
426               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
427               ! -----------------------------------             ===========================
428#if defined key_z_first
429               DO jj = 1, jpj
430                  DO ji = 1, jpim1
431                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
432#else
433               DO jk = 1, jpkm1
434                  DO jj = 1, jpj
435                     DO ji = 1, jpim1
436#endif
437                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
438                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
439                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
440                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
441                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
442                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
443                        !
444                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
445                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
446                        !
447                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
448                        vb(ji,jj,jk) = zvf
449                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
450                        vn(ji,jj,jk) = va(ji,jj,jk)
451                     END DO
452                  END DO
453               END DO
454#if defined key_z_first
455               DO jj = 1, jpj
456                  DO ji = 1, jpim1
457                     DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
458                        fse3u_b(ji,jj,jk) = ze3u_f(ji,jj,jk)                   ! e3u_b <-- filtered scale factor
459                        fse3v_b(ji,jj,jk) = ze3v_f(ji,jj,jk)
460                     END DO
461                  END DO
462               END DO
463#else
464               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
465               fse3v_b(:,:,:) = ze3v_f(:,:,:)
466#endif
467               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
468               CALL lbc_lnk( vb, 'V', -1. )
469            ENDIF
470            !
471         ENDIF
472         !
473      ENDIF
474
475#if defined ARPDBGSUM
476      WRITE(*,*) narea,': ARPDBG: dyn_nxt end SUM(un)=',SUM(un),'at step=',kt
477#endif
478
479      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
480         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
481      !
482      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
483      !
484   END SUBROUTINE dyn_nxt
485
486   !!=========================================================================
487END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.