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

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3865

Last change on this file since 3865 was 3865, checked in by acc, 11 years ago

Branch 2013/dev_r3858_NOC_ZTC, #863. Nearly complete port of 2011/dev_r2739_LOCEAN8_ZTC development branch into v3.5aplha base. Compiles and runs but currently unstable after 8 timesteps with ORCA2_LIM reference configuration.
  • Property svn:keywords set to Id
File size: 14.3 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 wrk_nemo        ! Memory Allocation
43   USE prtctl          ! Print control
44#if defined key_agrif
45   USE agrif_opa_interp
46#endif
47   USE timing          ! Timing
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC    dyn_nxt   ! routine called by step.F90
53
54   !! * Substitutions
55#  include "domzgr_substitute.h90"
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
58   !! $Id$
59   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61CONTAINS
62
63   SUBROUTINE dyn_nxt ( kt )
64      !!----------------------------------------------------------------------
65      !!                  ***  ROUTINE dyn_nxt  ***
66      !!                   
67      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
68      !!             condition on the after velocity, achieved the time stepping
69      !!             by applying the Asselin filter on now fields and swapping
70      !!             the fields.
71      !!
72      !! ** Method  : * After velocity is compute using a leap-frog scheme:
73      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
74      !!             Note that with flux form advection and variable volume layer
75      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
76      !!             velocity.
77      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
78      !!             the time stepping has already been done in dynspg module
79      !!
80      !!              * Apply lateral boundary conditions on after velocity
81      !!             at the local domain boundaries through lbc_lnk call,
82      !!             at the one-way open boundaries (lk_obc=T),
83      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
84      !!
85      !!              * Apply the time filter applied and swap of the dynamics
86      !!             arrays to start the next time step:
87      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
88      !!                (un,vn) = (ua,va).
89      !!             Note that with flux form advection and variable volume layer
90      !!             (lk_vvl=T), the time filter is applied on thickness weighted
91      !!             velocity.
92      !!
93      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
94      !!               un,vn   now horizontal velocity of next time-step
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      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
106      !!----------------------------------------------------------------------
107      !
108      IF( nn_timing == 1 )  CALL timing_start('dyn_nxt')
109      !
110      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
111      !
112      IF( kt == nit000 ) THEN
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
115         IF(lwp) WRITE(numout,*) '~~~~~~~'
116      ENDIF
117
118#if defined key_dynspg_flt
119      !
120      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
121      ! -------------
122
123      ! Update after velocity on domain lateral boundaries      (only local domain required)
124      ! --------------------------------------------------
125      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
126      CALL lbc_lnk( va, 'V', -1. ) 
127      !
128#else
129      ! Next velocity :   Leap-frog time stepping
130      ! -------------
131      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
132      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
133      !
134      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
135         DO jk = 1, jpkm1
136            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
137            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
138         END DO
139      ELSE                                            ! applied on thickness weighted velocity
140         DO jk = 1, jpkm1
141            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
142               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
143               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
144            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
145               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
146               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
147         END DO
148      ENDIF
149
150
151      ! Update after velocity on domain lateral boundaries
152      ! --------------------------------------------------     
153      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
154      CALL lbc_lnk( va, 'V', -1. ) 
155      !
156# if defined key_obc
157      !                                !* OBC open boundaries
158      IF( lk_obc ) CALL obc_dyn( kt )
159      !
160      IF( .NOT. lk_dynspg_flt ) THEN
161         ! Flather boundary condition : - Update sea surface height on each open boundary
162         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
163         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
164         !                              - Correct the barotropic velocities
165         IF( lk_obc ) CALL obc_dyn_bt( kt )
166         !
167!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
168         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
169         !
170         IF( lk_obc .AND. ln_vol_cst )   CALL obc_vol( kt )
171         !
172         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
173      ENDIF
174      !
175# elif defined key_bdy
176      !                                !* BDY open boundaries
177      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
178      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
179
180!!$   Do we need a call to bdy_vol here??
181      !
182# endif
183      !
184# if defined key_agrif
185      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
186# endif
187#endif
188
189      ! Time filter and swap of dynamics arrays
190      ! ------------------------------------------
191      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
192         DO jk = 1, jpkm1
193            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
194            vn(:,:,jk) = va(:,:,jk)
195         END DO
196      ELSE                                             !* Leap-Frog : Asselin filter and swap
197         !                                ! =============!
198         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
199            !                             ! =============!
200            DO jk = 1, jpkm1                             
201               DO jj = 1, jpj
202                  DO ji = 1, jpi   
203                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )
204                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )
205                     !
206                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
207                     vb(ji,jj,jk) = zvf
208                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
209                     vn(ji,jj,jk) = va(ji,jj,jk)
210                  END DO
211               END DO
212            END DO
213            !                             ! ================!
214         ELSE                             ! Variable volume !
215            !                             ! ================!
216            ! Before scale factor at t-points
217            ! (used as a now filtered scale factor until the swap)
218            ! ----------------------------------------------------
219            fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
220            ! Add volume filter correction: compatibility with tracer advection scheme
221            ! => time filter + conservation correction (only at the first level)
222            fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
223            !
224            IF( ln_dynadv_vec ) THEN
225               ! Before scale factor at (u/v)-points
226               ! -----------------------------------
227               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
228               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
229               ! Leap-Frog - Asselin filter and swap: applied on velocity
230               ! -----------------------------------
231               DO jk = 1, jpkm1
232                  DO jj = 1, jpj
233                     DO ji = 1, jpi
234                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
235                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
236                        !
237                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
238                        vb(ji,jj,jk) = zvf
239                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
240                        vn(ji,jj,jk) = va(ji,jj,jk)
241                     END DO
242                  END DO
243               END DO
244               !
245            ELSE
246               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
247               !------------------------------------------------
248               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
249               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
250               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
251               ! -----------------------------------             ===========================
252               DO jk = 1, jpkm1
253                  DO jj = 1, jpj
254                     DO ji = 1, jpi                   !                              ---------------------------
255                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
256                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
257                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
258                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
259                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
260                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
261                        !
262                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
263                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
264                        !
265                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
266                        vb(ji,jj,jk) = zvf
267                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
268                        vn(ji,jj,jk) = va(ji,jj,jk)
269                     END DO
270                  END DO
271               END DO
272               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
273               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
274            ENDIF
275            !
276         ENDIF
277         !
278      ENDIF
279
280      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
281         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
282      !
283      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
284      !
285      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
286      !
287   END SUBROUTINE dyn_nxt
288
289   !!=========================================================================
290END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.