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

source: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2800

Last change on this file since 2800 was 2800, checked in by davestorkey, 13 years ago
  1. Application of boundary conditions to barotropic and baroclinic velocities clearly separated.
  2. Option to input full velocities in boundary data (default expects barotropic and baroclinic velocities separately).
  3. Option to use initial conditions as boundary conditions coded.
  • Property svn:keywords set to Id
File size: 15.6 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 baroclinic velocities
33   USE obcvol          ! ocean open boundary condition (obc_vol routines)
34   USE in_out_manager  ! I/O manager
35   USE lbclnk          ! lateral boundary condition (or mpp link)
36   USE lib_mpp         ! MPP library
37   USE prtctl          ! Print control
38#if defined key_agrif
39   USE agrif_opa_interp
40#endif
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC    dyn_nxt   ! routine called by step.F90
46
47   !! * Substitutions
48#  include "domzgr_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE dyn_nxt ( kt )
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE dyn_nxt  ***
59      !!                   
60      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
61      !!             condition on the after velocity, achieved the time stepping
62      !!             by applying the Asselin filter on now fields and swapping
63      !!             the fields.
64      !!
65      !! ** Method  : * After velocity is compute using a leap-frog scheme:
66      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
67      !!             Note that with flux form advection and variable volume layer
68      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
69      !!             velocity.
70      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
71      !!             the time stepping has already been done in dynspg module
72      !!
73      !!              * Apply lateral boundary conditions on after velocity
74      !!             at the local domain boundaries through lbc_lnk call,
75      !!             at the one-way open boundaries (lk_obc=T),
76      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
77      !!
78      !!              * Apply the time filter applied and swap of the dynamics
79      !!             arrays to start the next time step:
80      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
81      !!                (un,vn) = (ua,va).
82      !!             Note that with flux form advection and variable volume layer
83      !!             (lk_vvl=T), the time filter is applied on thickness weighted
84      !!             velocity.
85      !!
86      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
87      !!               un,vn   now horizontal velocity of next time-step
88      !!----------------------------------------------------------------------
89      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
90      USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
91      USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
92      !
93      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
94      !
95      INTEGER  ::   ji, jj, jk   ! dummy loop indices
96#if ! defined key_dynspg_flt
97      REAL(wp) ::   z2dt         ! temporary scalar
98#endif
99      REAL(wp) ::   zue3a, zue3n, zue3b, zuf    ! local scalars
100      REAL(wp) ::   zve3a, zve3n, zve3b, zvf    !   -      -
101      REAL(wp) ::   zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1
102      !!----------------------------------------------------------------------
103
104      IF( wrk_in_use(2, 1,2,3) ) THEN
105         CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable')   ;   RETURN
106      ENDIF
107
108      IF( kt == nit000 ) THEN
109         IF(lwp) WRITE(numout,*)
110         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
111         IF(lwp) WRITE(numout,*) '~~~~~~~'
112      ENDIF
113
114#if defined key_dynspg_flt
115      !
116      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
117      ! -------------
118
119      ! Update after velocity on domain lateral boundaries      (only local domain required)
120      ! --------------------------------------------------
121      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
122      CALL lbc_lnk( va, 'V', -1. ) 
123      !
124#else
125      ! Next velocity :   Leap-frog time stepping
126      ! -------------
127      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
128      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
129      !
130      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
131         DO jk = 1, jpkm1
132            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
133            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
134         END DO
135      ELSE                                            ! applied on thickness weighted velocity
136         DO jk = 1, jpkm1
137            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
138               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
139               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
140            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
141               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
142               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
143         END DO
144      ENDIF
145
146
147      ! Update after velocity on domain lateral boundaries
148      ! --------------------------------------------------     
149      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
150      CALL lbc_lnk( va, 'V', -1. ) 
151      !
152# if defined key_obc
153      !                                !* OBC open boundaries
154      IF( .NOT. lk_dynspg_flt ) THEN
155
156         CALL obc_dyn( kt )
157
158!!$!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
159!!$         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
160!!$         !
161!!$         IF( ln_vol_cst )   CALL obc_vol( kt )
162!!$         !
163!!$         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
164
165      ENDIF
166      !
167# endif
168      !
169# if defined key_agrif
170      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
171# endif
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      ELSE                                             !* Leap-Frog : Asselin filter and swap
182         !                                ! =============!
183         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
184            !                             ! =============!
185            DO jk = 1, jpkm1                             
186               DO jj = 1, jpj
187                  DO ji = 1, jpi   
188                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
189                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
190                     !
191                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
192                     vb(ji,jj,jk) = zvf
193                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
194                     vn(ji,jj,jk) = va(ji,jj,jk)
195                  END DO
196               END DO
197            END DO
198            !                             ! ================!
199         ELSE                             ! Variable volume !
200            !                             ! ================!
201            ! Before scale factor at t-points
202            ! -------------------------------
203            DO jk = 1, jpkm1
204               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
205                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
206                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
207            ENDDO
208            ! Add volume filter correction only at the first level of t-point scale factors
209            zec = atfp * rdt / rau0
210            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
211            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
212            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
213            zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
214            zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
215            !
216            IF( ln_dynadv_vec ) THEN
217               ! Before scale factor at (u/v)-points
218               ! -----------------------------------
219               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
220               DO jk = 1, jpkm1
221                  DO jj = 1, jpjm1
222                     DO ji = 1, jpim1
223                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
224                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
225                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
226                        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) )
227                        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) )
228                     END DO
229                  END DO
230               END DO
231               ! lateral boundary conditions
232               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
233               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
234               ! Add initial scale factor to scale factor anomaly
235               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
236               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
237               ! Leap-Frog - Asselin filter and swap: applied on velocity
238               ! -----------------------------------
239               DO jk = 1, jpkm1
240                  DO jj = 1, jpj
241                     DO ji = 1, jpi
242                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
243                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
244                        !
245                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
246                        vb(ji,jj,jk) = zvf
247                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
248                        vn(ji,jj,jk) = va(ji,jj,jk)
249                     END DO
250                  END DO
251               END DO
252               !
253            ELSE
254               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
255               !-----------------------------------------------
256               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
257               DO jk = 1, jpkm1
258                  DO jj = 1, jpjm1
259                     DO ji = 1, jpim1
260                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
261                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
262                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
263                        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) )
264                        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) )
265                     END DO
266                  END DO
267               END DO
268               ! lateral boundary conditions
269               CALL lbc_lnk( ze3u_f, 'U', 1. )
270               CALL lbc_lnk( ze3v_f, 'V', 1. )
271               ! Add initial scale factor to scale factor anomaly
272               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
273               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
274               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
275               ! -----------------------------------             ===========================
276               DO jk = 1, jpkm1
277                  DO jj = 1, jpj
278                     DO ji = 1, jpim1
279                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
280                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
281                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
282                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
283                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
284                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
285                        !
286                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
287                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
288                        !
289                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
290                        vb(ji,jj,jk) = zvf
291                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
292                        vn(ji,jj,jk) = va(ji,jj,jk)
293                     END DO
294                  END DO
295               END DO
296               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
297               fse3v_b(:,:,:) = ze3v_f(:,:,:)
298               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
299               CALL lbc_lnk( vb, 'V', -1. )
300            ENDIF
301            !
302         ENDIF
303         !
304      ENDIF
305
306      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
307         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
308      !
309      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
310      !
311   END SUBROUTINE dyn_nxt
312
313   !!=========================================================================
314END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.