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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3317

Last change on this file since 3317 was 3317, checked in by gm, 12 years ago

Ediag branche: #927 restructuration of the trdicp computation - part I

  • Property svn:keywords set to Id
File size: 15.6 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   !!            3.4  !  2012-02  (G. Madec) add the diagnostic of the time filter trends
20   !!-------------------------------------------------------------------------
21 
22   !!-------------------------------------------------------------------------
23   !!   dyn_nxt      : obtain the next (after) horizontal velocity
24   !!-------------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE sbc_oce         ! Surface boundary condition: ocean fields
28   USE trdmod_oce      ! ocean trends
29   USE phycst          ! physical constants
30   USE dynspg_oce      ! type of surface pressure gradient
31   USE dynadv          ! dynamics: vector invariant versus flux form
32   USE domvvl          ! variable volume
33   USE trdmod         ! ocean dynamics trends
34   USE trdmod_oce     ! ocean variables trends
35   USE obc_oce         ! ocean open boundary conditions
36   USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
37   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
38   USE obcvol          ! ocean open boundary condition (obc_vol routines)
39   USE bdy_oce         ! ocean open boundary conditions
40   USE bdydta          ! ocean open boundary conditions
41   USE bdydyn          ! ocean open boundary conditions
42   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
43   USE in_out_manager  ! I/O manager
44   USE iom             ! I/O manager library
45   USE lbclnk          ! lateral boundary condition (or mpp link)
46   USE lib_mpp         ! MPP library
47   USE wrk_nemo        ! Memory Allocation
48   USE prtctl          ! Print control
49   USE timing          ! Timing
50#if defined key_agrif
51   USE agrif_opa_interp
52#endif
53
54   IMPLICIT NONE
55   PRIVATE
56
57   PUBLIC    dyn_nxt   ! routine called by step.F90
58
59   !! * Substitutions
60#  include "domzgr_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
63   !! $Id$
64   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   SUBROUTINE dyn_nxt ( kt )
69      !!----------------------------------------------------------------------
70      !!                  ***  ROUTINE dyn_nxt  ***
71      !!                   
72      !! ** Purpose :   Compute the after horizontal velocity. Apply the boundary
73      !!             condition on the after velocity, achieved the time stepping
74      !!             by applying the Asselin filter on now fields and swapping
75      !!             the fields.
76      !!
77      !! ** Method  : * After velocity is compute using a leap-frog scheme:
78      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
79      !!             Note that with flux form advection and variable volume layer
80      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
81      !!             velocity.
82      !!             Note also that in filtered free surface (lk_dynspg_flt=T),
83      !!             the time stepping has already been done in dynspg module
84      !!
85      !!              * Apply lateral boundary conditions on after velocity
86      !!             at the local domain boundaries through lbc_lnk call,
87      !!             at the one-way open boundaries (lk_obc=T),
88      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
89      !!
90      !!              * Apply the time filter applied and swap of the dynamics
91      !!             arrays to start the next time step:
92      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
93      !!                (un,vn) = (ua,va).
94      !!             Note that with flux form advection and variable volume layer
95      !!             (lk_vvl=T), the time filter is applied on thickness weighted
96      !!             velocity.
97      !!
98      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
99      !!               un,vn   now horizontal velocity of next time-step
100      !!----------------------------------------------------------------------
101      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
102      !
103      INTEGER  ::   ji, jj, jk   ! dummy loop indices
104      INTEGER  ::   iku, ikv     ! local integers
105#if ! defined key_dynspg_flt
106      REAL(wp) ::   z2dt         ! temporary scalar
107#endif
108      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars
109      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
110      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
111      !!----------------------------------------------------------------------
112      !
113      IF( nn_timing == 1 )  CALL timing_start('dyn_nxt')
114      !
115      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva )
116      !
117      IF( kt == nit000 ) THEN
118         IF(lwp) WRITE(numout,*)
119         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
120         IF(lwp) WRITE(numout,*) '~~~~~~~'
121      ENDIF
122
123#if defined key_dynspg_flt
124      !
125      ! Next velocity :   Leap-frog time stepping already done in dynspg_flt.F routine
126      ! -------------
127
128      ! Update after velocity on domain lateral boundaries      (only local domain required)
129      ! --------------------------------------------------
130      CALL lbc_lnk( ua, 'U', -1. )         ! local domain boundaries
131      CALL lbc_lnk( va, 'V', -1. ) 
132      !
133#else
134      ! Next velocity :   Leap-frog time stepping
135      ! -------------
136      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
137      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
138      !
139      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
140         DO jk = 1, jpkm1
141            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
142            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
143         END DO
144      ELSE                                            ! applied on thickness weighted velocity
145         DO jk = 1, jpkm1
146            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
147               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
148               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
149            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
150               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
151               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
152         END DO
153      ENDIF
154
155
156      ! Update after velocity on domain lateral boundaries
157      ! --------------------------------------------------     
158      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
159      CALL lbc_lnk( va, 'V', -1. ) 
160      !
161# if defined key_obc
162      !                                !* OBC open boundaries
163      CALL obc_dyn( kt )
164      !
165      IF( .NOT. lk_dynspg_flt ) THEN
166         ! Flather boundary condition : - Update sea surface height on each open boundary
167         !                                       sshn   (= after ssh   ) for explicit case (lk_dynspg_exp=T)
168         !                                       sshn_b (= after ssha_b) for time-splitting case (lk_dynspg_ts=T)
169         !                              - Correct the barotropic velocities
170         CALL obc_dyn_bt( kt )
171         !
172!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
173         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
174         !
175         IF( ln_vol_cst )   CALL obc_vol( kt )
176         !
177         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
178      ENDIF
179      !
180# elif defined key_bdy
181      !                                !* BDY open boundaries
182      IF( lk_dynspg_exp ) CALL bdy_dyn( kt )
183      IF( lk_dynspg_ts )  CALL bdy_dyn( kt, dyn3d_only=.true. )
184
185!!$   Do we need a call to bdy_vol here??
186      !
187# endif
188      !
189# if defined key_agrif
190      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
191# endif
192#endif
193
194      IF( ln_3D_trd_d ) THEN             ! 3D output: total momentum trends a prepare the atf trend computation
195         z1_2dt = 1._wp / (2. * rdt)            ! Euler or leap-frog time step
196         IF( neuler == 0 .AND. kt == nit000 )  z1_2dt = 1._wp / rdt
197         zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
198         zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
199         CALL iom_put( "utrd_tot", zua )        ! total momentum trends (but the asselin time filter)
200         CALL iom_put( "vtrd_tot", zva )
201         zua(:,:,:) = un(:,:,:)                 ! save the before velocity before the asselin filter
202         zva(:,:,:) = vn(:,:,:)                 ! (caution: there is a shift by 1 timestep in the
203         !                                      !  computation of the asselin filter trends)
204      ENDIF
205
206      ! Time filter and swap of dynamics arrays
207      ! ------------------------------------------
208      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
209         DO jk = 1, jpkm1
210            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
211            vn(:,:,jk) = va(:,:,jk)
212         END DO
213      ELSE                                             !* Leap-Frog : Asselin filter and swap
214         !                                ! =============!
215         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
216            !                             ! =============!
217            DO jk = 1, jpkm1                             
218               DO jj = 1, jpj
219                  DO ji = 1, jpi   
220                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
221                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
222                     !
223                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
224                     vb(ji,jj,jk) = zvf
225                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
226                     vn(ji,jj,jk) = va(ji,jj,jk)
227                  END DO
228               END DO
229            END DO
230            !                             ! ================!
231         ELSE                             ! Variable volume !
232            !                             ! ================!
233            !
234            DO jk = 1, jpkm1                 ! Before scale factor at t-points
235               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)                                   &
236                  &              + atfp * (  fse3t_b(:,:,jk) + fse3t_a(:,:,jk)     &
237                  &                         - 2._wp * fse3t_n(:,:,jk)            )
238            END DO
239            zec = atfp * rdt / rau0          ! Add filter correction only at the 1st level of t-point scale factors
240            fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
241            !
242            IF( ln_dynadv_vec ) THEN         ! vector invariant form (no thickness weighted calulation)
243               !
244               !                                      ! before scale factors at u- & v-pts (computed from fse3t_b)
245               CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) )
246               !
247               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap: applied on velocity
248                  DO jj = 1, jpj                      !                                                 --------
249                     DO ji = 1, jpi
250                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )
251                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )
252                        !
253                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
254                        vb(ji,jj,jk) = zvf
255                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
256                        vn(ji,jj,jk) = va(ji,jj,jk)
257                     END DO
258                  END DO
259               END DO
260               !
261            ELSE                             ! flux form (thickness weighted calulation)
262               !
263               CALL dom_vvl_2( kt, ze3u_f, ze3v_f )   ! before scale factors at u- & v-pts (computed from fse3t_b)
264               !
265               DO jk = 1, jpkm1                       ! Leap-Frog - Asselin filter and swap:
266                  DO jj = 1, jpj                      !                   applied on thickness weighted velocity
267                     DO ji = 1, jpim1                 !                              ---------------------------
268                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
269                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
270                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
271                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
272                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
273                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
274                        !
275                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
276                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
277                        !
278                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
279                        vb(ji,jj,jk) = zvf
280                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
281                        vn(ji,jj,jk) = va(ji,jj,jk)
282                     END DO
283                  END DO
284               END DO
285               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
286               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
287               CALL lbc_lnk( ub, 'U', -1. )                    ! lateral boundary conditions
288               CALL lbc_lnk( vb, 'V', -1. )
289            ENDIF
290            !
291         ENDIF
292         !
293      ENDIF
294
295      IF( ln_3D_trd_d ) THEN             ! 3D output: asselin filter trends on momentum
296         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
297         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
298         CALL trd_mod( zua, zva, jpdyn_trd_atf, 'DYN', kt )
299      ENDIF
300
301      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
302         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
303      !
304      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva )
305      !
306      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
307      !
308   END SUBROUTINE dyn_nxt
309
310   !!=========================================================================
311END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.