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

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