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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2723

Last change on this file since 2723 was 2723, checked in by rblod, 13 years ago

Correction in dynnxt for bdy, see ticket #811
  • Property svn:keywords set to Id
File size: 17.5 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   !! * Substitutions
52#  include "domzgr_substitute.h90"
53   !!----------------------------------------------------------------------
54   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
55   !! $Id$
56   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
57   !!----------------------------------------------------------------------
58CONTAINS
59
60   SUBROUTINE dyn_nxt ( kt )
61      !!----------------------------------------------------------------------
62      !!                  ***  ROUTINE dyn_nxt  ***
63      !!                   
64      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
65      !!             condition on the after velocity, achieved the time stepping
66      !!             by applying the Asselin filter on now fields and swapping
67      !!             the fields.
68      !!
69      !! ** Method  : * After velocity is compute using a leap-frog scheme:
70      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
71      !!             Note that with flux form advection and variable volume layer
72      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
73      !!             velocity.
74      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
75      !!             the time stepping has already been done in dynspg module
76      !!
77      !!              * Apply lateral boundary conditions on after velocity
78      !!             at the local domain boundaries through lbc_lnk call,
79      !!             at the radiative open boundaries (lk_obc=T),
80      !!             at the relaxed   open boundaries (lk_bdy=T), and
81      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
82      !!
83      !!              * Apply the time filter applied and swap of the dynamics
84      !!             arrays to start the next time step:
85      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
86      !!                (un,vn) = (ua,va).
87      !!             Note that with flux form advection and variable volume layer
88      !!             (lk_vvl=T), the time filter is applied on thickness weighted
89      !!             velocity.
90      !!
91      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
92      !!               un,vn   now horizontal velocity of next time-step
93      !!----------------------------------------------------------------------
94      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
95      USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
96      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
97      !
98      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
99      !
100      INTEGER  ::   ji, jj, jk   ! dummy loop indices
101#if ! defined key_dynspg_flt
102      REAL(wp) ::   z2dt         ! temporary scalar
103#endif
104      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
105      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
106      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
107      !!----------------------------------------------------------------------
108
109      IF( wrk_in_use(2, 1,2,3) ) THEN
110         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
111      ENDIF
112
113      IF( kt == nit000 ) THEN
114         IF(lwp) WRITE(numout,*)
115         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
116         IF(lwp) WRITE(numout,*) '~~~~~~~'
117      ENDIF
118
119#if defined key_dynspg_flt
120      !
121      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
122      ! -------------
123
124      ! Update after velocity on domain lateral boundaries      (only local domain required)
125      ! --------------------------------------------------
126      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
127      CALL lbc_lnk( va, 'V', -1. ) 
128      !
129#else
130      ! Next velocity :   Leap-frog time stepping
131      ! -------------
132      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
133      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
134      !
135      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
136         DO jk = 1, jpkm1
137            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
138            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
139         END DO
140      ELSE                                            ! applied on thickness weighted velocity
141         DO jk = 1, jpkm1
142            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
143               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
144               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
145            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
146               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
147               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
148         END DO
149      ENDIF
150
151
152      ! Update after velocity on domain lateral boundaries
153      ! --------------------------------------------------     
154      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
155      CALL lbc_lnk( va, 'V', -1. ) 
156      !
157# if defined key_obc
158      !                                !* OBC open boundaries
159      CALL obc_dyn( kt )
160      !
161      IF( .NOT. lk_dynspg_flt ) THEN
162         ! Flather boundary condition : - Update sea surface height on each open boundary
163         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
164         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
165         !                              - Correct the barotropic velocities
166         CALL obc_dyn_bt( kt )
167         !
168!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
169         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
170         !
171         IF( ln_vol_cst )   CALL obc_vol( kt )
172         !
173         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
174      ENDIF
175      !
176# elif defined key_bdy 
177      !                                !* BDY open boundaries
178      IF( .NOT. lk_dynspg_flt ) THEN
179         CALL bdy_dyn_frs( kt )
180#  if ! defined key_vvl
181         ua_e(:,:) = 0.e0
182         va_e(:,:) = 0.e0
183         ! Set these variables for use in bdy_dyn_fla
184         hur_e(:,:) = hur(:,:)
185         hvr_e(:,:) = hvr(:,:)
186         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
187            ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
188            va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
189         END DO
190         ua_e(:,:) = ua_e(:,:) * hur(:,:)
191         va_e(:,:) = va_e(:,:) * hvr(:,:)
192         DO jk = 1 , jpkm1
193            ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
194            va(:,:,jk) = va(:,:,jk) - va_e(:,:)
195         END DO
196         CALL bdy_dta_fla( kt+1, 0,2*nn_baro)
197         CALL bdy_dyn_fla( sshn_b )
198         CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
199         CALL lbc_lnk( va_e, 'V', -1. )   !
200         DO jk = 1 , jpkm1
201            ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
202            va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
203         END DO
204#  endif
205      ENDIF
206# endif
207      !
208# if defined key_agrif
209      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
210# endif
211#endif
212
213      ! Time filter and swap of dynamics arrays
214      ! ------------------------------------------
215      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
216         DO jk = 1, jpkm1
217            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
218            vn(:,:,jk) = va(:,:,jk)
219         END DO
220      ELSE                                             !* Leap-Frog : Asselin filter and swap
221         !                                ! =============!
222         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
223            !                             ! =============!
224            DO jk = 1, jpkm1                             
225               DO jj = 1, jpj
226                  DO ji = 1, jpi   
227                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
228                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
229                     !
230                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
231                     vb(ji,jj,jk) = zvf
232                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
233                     vn(ji,jj,jk) = va(ji,jj,jk)
234                  END DO
235               END DO
236            END DO
237            !                             ! ================!
238         ELSE                             ! Variable volume !
239            !                             ! ================!
240            ! Before scale factor at t-points
241            ! -------------------------------
242            DO jk = 1, jpkm1
243               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
244                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
245                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
246            ENDDO
247            ! Add volume filter correction only at the first level of t-point scale factors
248            zec = atfp * rdt / rau0
249            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
250            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
251            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
252            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
253            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
254            !
255            IF( ln_dynadv_vec ) THEN
256               ! Before scale factor at (u/v)-points
257               ! -----------------------------------
258               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
259               DO jk = 1, jpkm1
260                  DO jj = 1, jpjm1
261                     DO ji = 1, jpim1
262                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
263                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
264                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
265                        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) )
266                        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) )
267                     END DO
268                  END DO
269               END DO
270               ! lateral boundary conditions
271               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
272               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
273               ! Add initial scale factor to scale factor anomaly
274               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
275               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
276               ! Leap-Frog - Asselin filter and swap: applied on velocity
277               ! -----------------------------------
278               DO jk = 1, jpkm1
279                  DO jj = 1, jpj
280                     DO ji = 1, jpi
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
293               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
294               !-----------------------------------------------
295               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
296               DO jk = 1, jpkm1
297                  DO jj = 1, jpjm1
298                     DO ji = 1, jpim1
299                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
300                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
301                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
302                        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) )
303                        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) )
304                     END DO
305                  END DO
306               END DO
307               ! lateral boundary conditions
308               CALL lbc_lnk( ze3u_f, 'U', 1. )
309               CALL lbc_lnk( ze3v_f, 'V', 1. )
310               ! Add initial scale factor to scale factor anomaly
311               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
312               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
313               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
314               ! -----------------------------------             ===========================
315               DO jk = 1, jpkm1
316                  DO jj = 1, jpj
317                     DO ji = 1, jpim1
318                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
319                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
320                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
321                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
322                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
323                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
324                        !
325                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
326                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
327                        !
328                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
329                        vb(ji,jj,jk) = zvf
330                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
331                        vn(ji,jj,jk) = va(ji,jj,jk)
332                     END DO
333                  END DO
334               END DO
335               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
336               fse3v_b(:,:,:) = ze3v_f(:,:,:)
337               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
338               CALL lbc_lnk( vb, 'V', -1. )
339            ENDIF
340            !
341         ENDIF
342         !
343      ENDIF
344
345      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
346         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
347      !
348      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
349      !
350   END SUBROUTINE dyn_nxt
351
352   !!=========================================================================
353END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.