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

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

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