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

source: branches/2011/dev_UKM0_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3062

Last change on this file since 3062 was 3062, checked in by rfurner, 12 years ago

ticket #885. added in changes from branches/2011/UKMO_MERCATOR_obc_bdy_merge@2888
  • Property svn:keywords set to Id
File size: 16.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 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         ! ocean open boundary conditions
36   USE bdydta          ! ocean open boundary conditions
37   USE bdydyn          ! ocean open boundary conditions
38   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
39   USE in_out_manager  ! I/O manager
40   USE lbclnk          ! lateral boundary condition (or mpp link)
41   USE lib_mpp         ! MPP library
42   USE prtctl          ! Print control
43#if defined key_agrif
44   USE agrif_opa_interp
45#endif
46
47   IMPLICIT NONE
48   PRIVATE
49
50   PUBLIC    dyn_nxt   ! routine called by step.F90
51
52   !! * Substitutions
53#  include "domzgr_substitute.h90"
54   !!----------------------------------------------------------------------
55   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
56   !! $Id$
57   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE dyn_nxt ( kt )
62      !!----------------------------------------------------------------------
63      !!                  ***  ROUTINE dyn_nxt  ***
64      !!                   
65      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
66      !!             condition on the after velocity, achieved the time stepping
67      !!             by applying the Asselin filter on now fields and swapping
68      !!             the fields.
69      !!
70      !! ** Method  : * After velocity is compute using a leap-frog scheme:
71      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
72      !!             Note that with flux form advection and variable volume layer
73      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
74      !!             velocity.
75      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
76      !!             the time stepping has already been done in dynspg module
77      !!
78      !!              * Apply lateral boundary conditions on after velocity
79      !!             at the local domain boundaries through lbc_lnk call,
80      !!             at the one-way open boundaries (lk_obc=T),
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( lk_dynspg_exp ) CALL bdy_dyn( kt )
179      IF( lk_dynspg_ts )  CALL bdy_dyn( kt, dyn3d_only=.true. )
180
181!!$   Do we need a call to bdy_vol here??
182      !
183# endif
184      !
185# if defined key_agrif
186      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
187# endif
188#endif
189
190      ! Time filter and swap of dynamics arrays
191      ! ------------------------------------------
192      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
193         DO jk = 1, jpkm1
194            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
195            vn(:,:,jk) = va(:,:,jk)
196         END DO
197      ELSE                                             !* Leap-Frog : Asselin filter and swap
198         !                                ! =============!
199         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
200            !                             ! =============!
201            DO jk = 1, jpkm1                             
202               DO jj = 1, jpj
203                  DO ji = 1, jpi   
204                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
205                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
206                     !
207                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
208                     vb(ji,jj,jk) = zvf
209                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
210                     vn(ji,jj,jk) = va(ji,jj,jk)
211                  END DO
212               END DO
213            END DO
214            !                             ! ================!
215         ELSE                             ! Variable volume !
216            !                             ! ================!
217            ! Before scale factor at t-points
218            ! -------------------------------
219            DO jk = 1, jpkm1
220               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
221                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
222                  &                         - 2.e0 * fse3t_n(:,:,jk)            )
223            ENDDO
224            ! Add volume filter correction only at the first level of t-point scale factors
225            zec = atfp * rdt / rau0
226            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
227            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
228            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
229            zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )
230            zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )
231            !
232            IF( ln_dynadv_vec ) THEN
233               ! Before scale factor at (u/v)-points
234               ! -----------------------------------
235               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
236               DO jk = 1, jpkm1
237                  DO jj = 1, jpjm1
238                     DO ji = 1, jpim1
239                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
240                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
241                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
242                        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) )
243                        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) )
244                     END DO
245                  END DO
246               END DO
247               ! lateral boundary conditions
248               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
249               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
250               ! Add initial scale factor to scale factor anomaly
251               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
252               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
253               ! Leap-Frog - Asselin filter and swap: applied on velocity
254               ! -----------------------------------
255               DO jk = 1, jpkm1
256                  DO jj = 1, jpj
257                     DO ji = 1, jpi
258                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
259                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
260                        !
261                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
262                        vb(ji,jj,jk) = zvf
263                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
264                        vn(ji,jj,jk) = va(ji,jj,jk)
265                     END DO
266                  END DO
267               END DO
268               !
269            ELSE
270               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
271               !-----------------------------------------------
272               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
273               DO jk = 1, jpkm1
274                  DO jj = 1, jpjm1
275                     DO ji = 1, jpim1
276                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
277                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
278                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
279                        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) )
280                        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) )
281                     END DO
282                  END DO
283               END DO
284               ! lateral boundary conditions
285               CALL lbc_lnk( ze3u_f, 'U', 1. )
286               CALL lbc_lnk( ze3v_f, 'V', 1. )
287               ! Add initial scale factor to scale factor anomaly
288               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
289               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
290               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
291               ! -----------------------------------             ===========================
292               DO jk = 1, jpkm1
293                  DO jj = 1, jpj
294                     DO ji = 1, jpim1
295                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
296                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
297                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
298                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
299                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
300                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
301                        !
302                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
303                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
304                        !
305                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
306                        vb(ji,jj,jk) = zvf
307                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
308                        vn(ji,jj,jk) = va(ji,jj,jk)
309                     END DO
310                  END DO
311               END DO
312               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
313               fse3v_b(:,:,:) = ze3v_f(:,:,:)
314               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
315               CALL lbc_lnk( vb, 'V', -1. )
316            ENDIF
317            !
318         ENDIF
319         !
320      ENDIF
321
322      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
323         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
324      !
325      IF( wrk_not_released(2, 1,2,3) )   CALL ctl_stop('dyn_nxt: failed to release workspace arrays')
326      !
327   END SUBROUTINE dyn_nxt
328
329   !!=========================================================================
330END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.