source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 4 years ago

All wrk_alloc removed

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