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/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 6225

Last change on this file since 6225 was 6225, checked in by jamesharle, 8 years ago

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

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