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

source: branches/2011/dev_r2802_NOCS_vvlfix/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2807

Last change on this file since 2807 was 2807, checked in by acc, 13 years ago

Branch: dev_r2802_NOCS_vvlfix. Bugfix corrections to the calculation of fse3u_b and fse3v_b to address problems with vvl and partial steps. See comments added to ticket #812
  • Property svn:keywords set to Id
File size: 15.0 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 oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
95      !
96      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
97      !
98      INTEGER  ::   ji, jj, jk   ! dummy loop indices
99      INTEGER  ::   iku, ikv     ! local integers
100#if ! defined key_dynspg_flt
101      REAL(wp) ::   z2dt         ! temporary scalar
102#endif
103      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
104      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
105      !!----------------------------------------------------------------------
106
107      IF( kt == nit000 ) THEN
108         IF(lwp) WRITE(numout,*)
109         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
110         IF(lwp) WRITE(numout,*) '~~~~~~~'
111      ENDIF
112
113#if defined key_dynspg_flt
114      !
115      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
116      ! -------------
117
118      ! Update after velocity on domain lateral boundaries      (only local domain required)
119      ! --------------------------------------------------
120      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
121      CALL lbc_lnk( va, 'V', -1. ) 
122      !
123#else
124      ! Next velocity :   Leap-frog time stepping
125      ! -------------
126      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
127      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
128      !
129      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
130         DO jk = 1, jpkm1
131            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
132            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
133         END DO
134      ELSE                                            ! applied on thickness weighted velocity
135         DO jk = 1, jpkm1
136            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
137               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
138               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
139            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
140               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
141               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
142         END DO
143      ENDIF
144
145
146      ! Update after velocity on domain lateral boundaries
147      ! --------------------------------------------------     
148      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
149      CALL lbc_lnk( va, 'V', -1. ) 
150      !
151# if defined key_obc
152      !                                !* OBC open boundaries
153      CALL obc_dyn( kt )
154      !
155      IF( .NOT. lk_dynspg_flt ) THEN
156         ! Flather boundary condition : - Update sea surface height on each open boundary
157         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
158         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
159         !                              - Correct the barotropic velocities
160         CALL obc_dyn_bt( kt )
161         !
162!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
163         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
164         !
165         IF( ln_vol_cst )   CALL obc_vol( kt )
166         !
167         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
168      ENDIF
169      !
170# elif defined key_bdy 
171      !                                !* BDY open boundaries
172      IF( .NOT. lk_dynspg_flt ) THEN
173         CALL bdy_dyn_frs( kt )
174#  if ! defined key_vvl
175         ua_e(:,:) = 0.e0
176         va_e(:,:) = 0.e0
177         ! Set these variables for use in bdy_dyn_fla
178         hur_e(:,:) = hur(:,:)
179         hvr_e(:,:) = hvr(:,:)
180         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
181            ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
182            va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
183         END DO
184         ua_e(:,:) = ua_e(:,:) * hur(:,:)
185         va_e(:,:) = va_e(:,:) * hvr(:,:)
186         DO jk = 1 , jpkm1
187            ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:)
188            va(:,:,jk) = va(:,:,jk) - va_e(:,:)
189         END DO
190         CALL bdy_dta_fla( kt+1, 0,2*nn_baro)
191         CALL bdy_dyn_fla( sshn_b )
192         CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
193         CALL lbc_lnk( va_e, 'V', -1. )   !
194         DO jk = 1 , jpkm1
195            ua(:,:,jk) = ( ua(:,:,jk) + ua_e(:,:) ) * umask(:,:,jk)
196            va(:,:,jk) = ( va(:,:,jk) + va_e(:,:) ) * vmask(:,:,jk)
197         END DO
198#  endif
199      ENDIF
200# endif
201      !
202# if defined key_agrif
203      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
204# endif
205#endif
206
207      ! Time filter and swap of dynamics arrays
208      ! ------------------------------------------
209      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
210         DO jk = 1, jpkm1
211            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
212            vn(:,:,jk) = va(:,:,jk)
213         END DO
214      ELSE                                             !* Leap-Frog : Asselin filter and swap
215         !                                ! =============!
216         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
217            !                             ! =============!
218            DO jk = 1, jpkm1                             
219               DO jj = 1, jpj
220                  DO ji = 1, jpi   
221                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
222                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
223                     !
224                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
225                     vb(ji,jj,jk) = zvf
226                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
227                     vn(ji,jj,jk) = va(ji,jj,jk)
228                  END DO
229               END DO
230            END DO
231            !                             ! ================!
232         ELSE                             ! Variable volume !
233            !                             ! ================!
234            !
235            DO jk = 1, jpkm1                 ! Before scale factor at t-points
236               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
237                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
238                  &                         - 2._wp * fse3t_n(:,:,jk)            )
239            END DO
240            zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors
241            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
242            !
243            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation)
244               !
245               !                                      ! before scale factors at u- & v-pts (computed from fse3t_b)
246               CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) )
247               !
248               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity
249                  DO jj = 1, jpj                      !                                                 --------
250                     DO ji = 1, jpi
251                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
252                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
253                        !
254                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
255                        vb(ji,jj,jk) = zvf
256                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
257                        vn(ji,jj,jk) = va(ji,jj,jk)
258                     END DO
259                  END DO
260               END DO
261               !
262            ELSE                             ! flux form (thickness weighted calulation)
263               !
264               CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b)
265               !
266               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap:
267                  DO jj = 1, jpj                      !                   applied on thickness weighted velocity
268                     DO ji = 1, jpim1                 !                              ---------------------------
269                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
270                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
271                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
272                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
273                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
274                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
275                        !
276                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
277                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
278                        !
279                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
280                        vb(ji,jj,jk) = zvf
281                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
282                        vn(ji,jj,jk) = va(ji,jj,jk)
283                     END DO
284                  END DO
285               END DO
286               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
287               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
288               CALL lbc_lnk( ub, 'U', -1. )                    ! lateral boundary conditions
289               CALL lbc_lnk( vb, 'V', -1. )
290            ENDIF
291            !
292         ENDIF
293         !
294      ENDIF
295
296      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
297         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
298      !
299   END SUBROUTINE dyn_nxt
300
301   !!=========================================================================
302END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.