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 @ 3318

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

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