source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 6079

Last change on this file since 6079 was 6079, checked in by jamesharle, 6 years ago

merge to trunk@5936

  • Property svn:keywords set to Id
File size: 17.9 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.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
20   !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
21   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification
22   !!-------------------------------------------------------------------------
23 
24   !!-------------------------------------------------------------------------
25   !!   dyn_nxt      : obtain the next (after) horizontal velocity
26   !!-------------------------------------------------------------------------
27   USE oce             ! ocean dynamics and tracers
28   USE dom_oce         ! ocean space and time domain
29   USE sbc_oce         ! Surface boundary condition: ocean fields
30   USE phycst          ! physical constants
31   USE dynadv          ! dynamics: vector invariant versus flux form
32   USE domvvl          ! variable volume
33   USE bdy_oce         ! ocean open boundary conditions
34   USE bdydta          ! ocean open boundary conditions
35   USE bdydyn          ! ocean open boundary conditions
36   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
37   USE trd_oce         ! trends: ocean variables
38   USE trddyn          ! trend manager: dynamics
39   USE trdken          ! trend manager: kinetic energy
40   !
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 :   Finalize after horizontal velocity. Apply the boundary
71      !!             condition on the after velocity, achieve the time stepping
72      !!             by applying the Asselin filter on now fields and swapping
73      !!             the fields.
74      !!
75      !! ** Method  : * Ensure after velocities transport matches time splitting
76      !!             estimate (ln_dynspg_ts=T)
77      !!
78      !!              * Apply lateral boundary conditions on after velocity
79      !!             at the local domain boundaries through lbc_lnk call,
80      !!             at the one-way open boundaries (lk_bdy=T),
81      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
82      !!
83      !!              * Apply the time filter applied and swap of the dynamics
84      !!             arrays to start the next time step:
85      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
86      !!                (un,vn) = (ua,va).
87      !!             Note that with flux form advection and variable volume layer
88      !!             (lk_vvl=T), the time filter is applied on thickness weighted
89      !!             velocity.
90      !!
91      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
92      !!               un,vn   now horizontal velocity of next time-step
93      !!----------------------------------------------------------------------
94      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
95      !
96      INTEGER  ::   ji, jj, jk   ! dummy loop indices
97      INTEGER  ::   iku, ikv     ! local integers
98      REAL(wp) ::   zue3a, zue3n, zue3b, zuf           ! local scalars
99      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
100      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
101      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
102      !!----------------------------------------------------------------------
103      !
104      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
105      !
106      IF( ln_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
107      IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
108      IF( l_trddyn     )   CALL wrk_alloc( jpi,jpj,jpk, zua, zva )
109      !
110      IF( kt == nit000 ) THEN
111         IF(lwp) WRITE(numout,*)
112         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
113         IF(lwp) WRITE(numout,*) '~~~~~~~'
114      ENDIF
115
116      IF ( ln_dynspg_ts ) THEN
117         ! Ensure below that barotropic velocities match time splitting estimate
118         ! Compute actual transport and replace it with ts estimate at "after" time step
119         zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
120         zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
121         DO jk = 2, jpkm1
122            zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
123            zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
124         END DO
125         DO jk = 1, jpkm1
126            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
127            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
128         END DO
129
130         IF ( .NOT.ln_bt_fw ) THEN
131            ! Remove advective velocity from "now velocities"
132            ! prior to asselin filtering     
133            ! In the forward case, this is done below after asselin filtering   
134            ! so that asselin contribution is removed at the same time
135            DO jk = 1, jpkm1
136               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
137               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
138            END DO 
139         ENDIF
140      ENDIF
141
142      ! Update after velocity on domain lateral boundaries
143      ! --------------------------------------------------     
144# if defined key_agrif
145      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
146# endif
147      !
148      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
149      CALL lbc_lnk( va, 'V', -1. ) 
150      !
151# if defined key_bdy
152      !                                !* BDY open boundaries
153      IF( lk_bdy .AND. ln_dynspg_exp ) CALL bdy_dyn( kt )
154      IF( lk_bdy .AND. ln_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
155
156!!$   Do we need a call to bdy_vol here??
157      !
158# endif
159      !
160      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
161         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
162         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
163         !
164         !                                  ! Kinetic energy and Conversion
165         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
166         !
167         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
168            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
169            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
170            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
171            CALL iom_put( "vtrd_tot", zva )
172         ENDIF
173         !
174         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
175         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
176         !                                  !  computation of the asselin filter trends)
177      ENDIF
178
179      ! Time filter and swap of dynamics arrays
180      ! ------------------------------------------
181      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
182         DO jk = 1, jpkm1
183            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
184            vn(:,:,jk) = va(:,:,jk)
185         END DO
186         IF (lk_vvl) THEN
187            DO jk = 1, jpkm1
188               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
189               fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
190               fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
191            ENDDO
192         ENDIF
193      ELSE                                             !* Leap-Frog : Asselin filter and swap
194         !                                ! =============!
195         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
196            !                             ! =============!
197            DO jk = 1, jpkm1                             
198               DO jj = 1, jpj
199                  DO ji = 1, jpi   
200                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
201                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
202                     !
203                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
204                     vb(ji,jj,jk) = zvf
205                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
206                     vn(ji,jj,jk) = va(ji,jj,jk)
207                  END DO
208               END DO
209            END DO
210            !                             ! ================!
211         ELSE                             ! Variable volume !
212            !                             ! ================!
213            ! Before scale factor at t-points
214            ! (used as a now filtered scale factor until the swap)
215            ! ----------------------------------------------------
216            IF (ln_dynspg_ts.AND.ln_bt_fw) THEN
217               ! No asselin filtering on thicknesses if forward time splitting
218               fse3t_b(:,:,:) = fse3t_n(:,:,:)
219            ELSE
220               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
221               ! Add volume filter correction: compatibility with tracer advection scheme
222               ! => time filter + conservation correction (only at the first level)
223               IF ( nn_isf == 0) THEN   ! if no ice shelf melting
224                  fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
225                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
226               ELSE                     ! if ice shelf melting
227                  DO jj = 1,jpj
228                     DO ji = 1,jpi
229                        jk = mikt(ji,jj)
230                        fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       &
231                                          &                          * ( (emp_b(ji,jj)    - emp(ji,jj)   ) &
232                                          &                            - (rnf_b(ji,jj)    - rnf(ji,jj)   ) &
233                                          &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk)
234                     END DO
235                  END DO
236               END IF
237            ENDIF
238            !
239            IF( ln_dynadv_vec ) THEN
240               ! Before scale factor at (u/v)-points
241               ! -----------------------------------
242               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
243               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
244               ! Leap-Frog - Asselin filter and swap: applied on velocity
245               ! -----------------------------------
246               DO jk = 1, jpkm1
247                  DO jj = 1, jpj
248                     DO ji = 1, jpi
249                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
250                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
251                        !
252                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
253                        vb(ji,jj,jk) = zvf
254                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
255                        vn(ji,jj,jk) = va(ji,jj,jk)
256                     END DO
257                  END DO
258               END DO
259               !
260            ELSE
261               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
262               !------------------------------------------------
263               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
264               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
265               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
266               ! -----------------------------------             ===========================
267               DO jk = 1, jpkm1
268                  DO jj = 1, jpj
269                     DO ji = 1, jpi                 
270                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
271                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
272                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
273                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
274                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
275                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
276                        !
277                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
278                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
279                        !
280                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
281                        vb(ji,jj,jk) = zvf
282                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
283                        vn(ji,jj,jk) = va(ji,jj,jk)
284                     END DO
285                  END DO
286               END DO
287               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
288               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
289            ENDIF
290            !
291         ENDIF
292         !
293         IF (ln_dynspg_ts.AND.ln_bt_fw) THEN
294            ! Revert "before" velocities to time split estimate
295            ! Doing it here also means that asselin filter contribution is removed 
296            zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
297            zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
298            DO jk = 2, jpkm1
299               zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
300               zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
301            END DO
302            DO jk = 1, jpkm1
303               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
304               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
305            END DO
306         ENDIF
307         !
308      ENDIF ! neuler =/0
309      !
310      ! Set "now" and "before" barotropic velocities for next time step:
311      ! JC: Would be more clever to swap variables than to make a full vertical
312      ! integration
313      !
314      !
315      IF (lk_vvl) THEN
316         hu_b(:,:) = 0.
317         hv_b(:,:) = 0.
318         DO jk = 1, jpkm1
319            hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
320            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
321         END DO
322         hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
323         hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
324      ENDIF
325      !
326      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
327      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
328      !
329      DO jk = 1, jpkm1
330         DO jj = 1, jpj
331            DO ji = 1, jpi
332               un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
333               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
334               !
335               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
336               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
337            END DO
338         END DO
339      END DO
340      !
341      !
342      un_b(:,:) = un_b(:,:) * hur_a(:,:)
343      vn_b(:,:) = vn_b(:,:) * hvr_a(:,:)
344      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
345      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
346      !
347      !
348
349      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
350         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
351         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
352         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
353      ENDIF
354      !
355      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
356         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
357      !
358      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
359      IF( lk_vvl.AND.(.NOT.ln_dynadv_vec ) ) CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
360      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk, zua, zva )
361      !
362      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
363      !
364   END SUBROUTINE dyn_nxt
365
366   !!=========================================================================
367END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.