source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 3876

Last change on this file since 3876 was 3876, checked in by gm, 8 years ago

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update

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