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

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3161

Last change on this file since 3161 was 3161, checked in by cetlod, 12 years ago

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