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

source: branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 2068

Last change on this file since 2068 was 2068, checked in by mlelod, 14 years ago

ticket: #663 ensuring restartability and conservation
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 17.4 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   !!-------------------------------------------------------------------------
18 
19   !!-------------------------------------------------------------------------
20   !!   dyn_nxt      : obtain the next (after) horizontal velocity
21   !!-------------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers
23   USE dom_oce         ! ocean space and time domain
24   USE sbc_oce         ! Surface boundary condition: ocean fields
25   USE phycst          ! physical constants
26   USE dynspg_oce      ! type of surface pressure gradient
27   USE dynadv          ! dynamics: vector invariant versus flux form
28   USE domvvl          ! variable volume
29   USE obc_oce         ! ocean open boundary conditions
30   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
31   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
32   USE obcvol          ! ocean open boundary condition (obc_vol routines)
33   USE bdy_oce         ! unstructured open boundary conditions
34   USE bdydta          ! unstructured open boundary conditions
35   USE bdydyn          ! unstructured open boundary conditions
36   USE agrif_opa_update
37   USE agrif_opa_interp
38   USE in_out_manager  ! I/O manager
39   USE lbclnk          ! lateral boundary condition (or mpp link)
40   USE prtctl          ! Print control
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.2 , LOCEAN-IPSL (2009)
51   !! $Id$
52   !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
53   !!-------------------------------------------------------------------------
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 radiative open boundaries (lk_obc=T),
77      !!             at the relaxed   open boundaries (lk_bdy=T), and
78      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
79      !!
80      !!              * Apply the time filter applied and swap of the dynamics
81      !!             arrays to start the next time step:
82      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
83      !!                (un,vn) = (ua,va).
84      !!             Note that with flux form advection and variable volume layer
85      !!             (lk_vvl=T), the time filter is applied on thickness weighted
86      !!             velocity.
87      !!
88      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
89      !!               un,vn   now horizontal velocity of next time-step
90      !!----------------------------------------------------------------------
91#if defined key_vvl
92      USE oce, ONLY :   ze3u_f => ta   ! use ta as 3D workspace
93      USE oce, ONLY :   ze3v_f => sa   ! use sa as 3D workspace
94#endif
95      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
96      !!
97      INTEGER  ::   ji, jj, jk   ! dummy loop indices
98#if ! defined key_dynspg_flt
99      REAL(wp) ::   z2dt         ! temporary scalar
100#endif
101      REAL(wp) ::   zue3a , zue3n , zue3b    ! temporary scalar
102      REAL(wp) ::   zve3a , zve3n , zve3b    !    -         -
103      REAL(wp) ::   zuf   , zvf              !    -         -
104      REAL(wp) ::   zec                      !    -         -
105      REAL(wp) ::   zv_t_ij  , zv_t_ip1j     !     -        -
106      REAL(wp) ::   zv_t_ijp1                !     -        -
107      REAL(wp), DIMENSION(jpi,jpj) ::  zs_t, zs_u_1, zs_v_1      ! temporary 2D workspace
108      !!----------------------------------------------------------------------
109
110      IF( kt == nit000 ) THEN
111         IF(lwp) WRITE(numout,*)
112         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
113         IF(lwp) WRITE(numout,*) '~~~~~~~'
114      ENDIF
115
116#if defined key_dynspg_flt
117      !
118      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
119      ! -------------
120
121      ! Update after velocity on domain lateral boundaries      (only local domain required)
122      ! --------------------------------------------------
123      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
124      CALL lbc_lnk( va, 'V', -1. ) 
125      !
126#else
127      ! Next velocity :   Leap-frog time stepping
128      ! -------------
129      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
130      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
131      !
132      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
133         DO jk = 1, jpkm1
134            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
135            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
136         END DO
137      ELSE                                            ! applied on thickness weighted velocity
138         DO jk = 1, jpkm1
139            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
140               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
141               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
142            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
143               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
144               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
145         END DO
146      ENDIF
147
148
149      ! Update after velocity on domain lateral boundaries
150      ! --------------------------------------------------     
151      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
152      CALL lbc_lnk( va, 'V', -1. ) 
153      !
154# if defined key_obc
155      !                                !* OBC open boundaries
156      CALL obc_dyn( kt )
157      !
158      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
159         ! Flather boundary condition : - Update sea surface height on each open boundary
160         !                                       sshn   (= after ssh   ) for explicit case
161         !                                       sshn_b (= after ssha_b) for time-splitting case
162         !                              - Correct the barotropic velocities
163         CALL obc_dyn_bt( kt )
164         !
165!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
166         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
167         !
168         IF( ln_vol_cst )   CALL obc_vol( kt )
169         !
170         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
171      ENDIF
172      !
173# elif defined key_bdy 
174      !                                !* BDY open boundaries
175      !RB all this part should be in a specific routine
176      IF( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN       ! except for filtered option
177         !
178         CALL bdy_dyn_frs( kt )
179         !
180         IF( ln_bdy_dyn_fla ) THEN
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_bt( kt+1, 0)
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         !
206      ENDIF
207# endif
208      !
209# if defined key_agrif
210      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
211# endif
212#endif
213
214      ! Time filter and swap of dynamics arrays
215      ! ------------------------------------------
216      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
217         DO jk = 1, jpkm1
218            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
219            vn(:,:,jk) = va(:,:,jk)
220         END DO
221      ELSE                                             !* Leap-Frog : Asselin filter and swap
222         !                                ! =============!
223         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
224            !                             ! =============!
225            DO jk = 1, jpkm1                             
226               DO jj = 1, jpj
227                  DO ji = 1, jpi   
228                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
229                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
230                     !
231                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
232                     vb(ji,jj,jk) = zvf
233                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
234                     vn(ji,jj,jk) = va(ji,jj,jk)
235                  END DO
236               END DO
237            END DO
238            !                             ! ================!
239         ELSE                             ! Variable volume !
240            !                             ! ================!
241            ! Before scale factor at t-points
242            ! -------------------------------
243            DO jk = 1, jpkm1
244               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) + atfp * fse3t_d(:,:,jk)
245            ENDDO
246            ! Add volume filter correction only at the first level of t-point scale factors
247            zec = atfp * rdt / rau0
248            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
249            ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations
250            zs_t  (:,:) =       e1t(:,:) * e2t(:,:)
251            zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:)
252            zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:)
253            !
254            IF( ln_dynadv_vec ) THEN
255               ! Before scale factor at (u/v)-points
256               ! -----------------------------------
257               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
258               DO jk = 1, jpkm1
259                  DO jj = 1, jpjm1
260                     DO ji = 1, jpim1
261                        zv_t_ij           = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
262                        zv_t_ip1j         = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
263                        zv_t_ijp1         = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
264                        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) )
265                        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) )
266                     END DO
267                  END DO
268               END DO
269               ! lateral boundary conditions
270               CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. )
271               CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. )
272               ! Add initial scale factor to scale factor anomaly
273               fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:)
274               fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:)
275               ! Leap-Frog - Asselin filter and swap: applied on velocity
276               ! -----------------------------------
277               DO jk = 1, jpkm1
278                  DO jj = 1, jpj
279                     DO ji = 1, jpi
280                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
281                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
282                        !
283                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
284                        vb(ji,jj,jk) = zvf
285                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
286                        vn(ji,jj,jk) = va(ji,jj,jk)
287                     END DO
288                  END DO
289               END DO
290               !
291            ELSE
292               ! Temporary filered scale factor at (u/v)-points (will become before scale factor)
293               !-----------------------------------------------
294               ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points
295               DO jk = 1, jpkm1
296                  DO jj = 1, jpjm1
297                     DO ji = 1, jpim1
298                        zv_t_ij          = zs_t(ji  ,jj  ) * fse3t_b(ji  ,jj  ,jk)
299                        zv_t_ip1j        = zs_t(ji+1,jj  ) * fse3t_b(ji+1,jj  ,jk)
300                        zv_t_ijp1        = zs_t(ji  ,jj+1) * fse3t_b(ji  ,jj+1,jk)
301                        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) )
302                        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) )
303                     END DO
304                  END DO
305               END DO
306               ! lateral boundary conditions
307               CALL lbc_lnk( ze3u_f, 'U', 1. )
308               CALL lbc_lnk( ze3v_f, 'V', 1. )
309               ! Add initial scale factor to scale factor anomaly
310               ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:)
311               ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:)
312               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
313               ! -----------------------------------             ===========================
314               DO jk = 1, jpkm1
315                  DO jj = 1, jpj
316                     DO ji = 1, jpim1
317                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
318                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
319                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
320                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
321                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
322                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
323                        !
324                        zuf  = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
325                        zvf  = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
326                        !
327                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
328                        vb(ji,jj,jk) = zvf
329                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
330                        vn(ji,jj,jk) = va(ji,jj,jk)
331                     END DO
332                  END DO
333               END DO
334               fse3u_b(:,:,:) = ze3u_f(:,:,:)                   ! e3u_b <-- filtered scale factor
335               fse3v_b(:,:,:) = ze3v_f(:,:,:)
336               CALL lbc_lnk( ub, 'U', -1. )                     ! lateral boundary conditions
337               CALL lbc_lnk( vb, 'V', -1. )
338            ENDIF
339            !
340         ENDIF
341         !
342      ENDIF
343
344#if defined key_agrif
345      ! Update velocity at AGRIF zoom boundaries
346      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Dyn( kt )
347#endif     
348
349      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
350         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
351      !
352   END SUBROUTINE dyn_nxt
353
354   !!=========================================================================
355END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.