source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynnxt.F90 @ 9939

Last change on this file since 9939 was 9939, checked in by gm, 2 years ago

#1911 (ENHANCE-04): RK3 branche phased with MLF@9937 branche

  • Property svn:keywords set to Id
File size: 18.2 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 sbcrnf         ! river runoffs
31   USE sbcisf         ! ice shelf
32   USE phycst         ! physical constants
33   USE dynadv         ! dynamics: vector invariant versus flux form
34   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
35   USE domvvl         ! variable volume
36   USE bdy_oce   , ONLY: ln_bdy
37   USE bdydta         ! ocean open boundary conditions
38   USE bdydyn         ! ocean open boundary conditions
39   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
40   USE trd_oce        ! trends: ocean variables
41   USE trddyn         ! trend manager: dynamics
42   USE trdken         ! trend manager: kinetic energy
43   !
44   USE in_out_manager ! I/O manager
45   USE iom            ! I/O manager library
46   USE lbclnk         ! lateral boundary condition (or mpp link)
47   USE lib_mpp        ! MPP library
48   USE prtctl         ! Print control
49   USE timing         ! Timing
50#if defined key_agrif
51   USE agrif_oce_interp
52#endif
53
54   IMPLICIT NONE
55   PRIVATE
56
57   PUBLIC    dyn_nxt   ! routine called by step.F90
58
59   !!----------------------------------------------------------------------
60   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
61   !! $Id$
62   !! Software governed by the CeCILL licence     (./LICENSE)
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 (ln_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) + rn_atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
86      !!                (un,vn) = (ua,va).
87      !!             Note that with flux form advection and non linear free surface,
88      !!             the time filter is applied on thickness weighted velocity.
89      !!             As a result, dyn_nxt MUST be called after tra_nxt.
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  ::   ikt          ! local integers
98      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef   ! local scalars
99      REAL(wp) ::   zve3a, zve3n, zve3b, zvf          !   -      -
100      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
101      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
102      !!----------------------------------------------------------------------
103      !
104      IF( ln_timing    )   CALL timing_start('dyn_nxt')
105      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
106      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
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(:,:)*r1_hu_n(:,:) + un_b(:,:) ) * umask(:,:,jk)
135               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + 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_multi( ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
147      !
148      !                                !* BDY open boundaries
149      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
150      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
151
152!!$   Do we need a call to bdy_vol here??
153      !
154      IF( l_trddyn ) THEN             !* prepare the atf trend computation + some diagnostics
155         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )    ! Kinetic energy and Conversion
156         IF( ln_dyn_trd ) THEN                                       ! 3D output: total momentum trends
157            IF( ln_dynadv_vec ) THEN
158               zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt
159               zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt
160            ELSE
161               zua(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt
162               zva(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt
163            ENDIF
164            CALL iom_put( "utrd_tot", zua )                          ! total momentum trends, except the asselin filter
165            CALL iom_put( "vtrd_tot", zva )
166         ENDIF
167         zua(:,:,:) = un(:,:,:)                    ! save the now velocity before the asselin filter
168         zva(:,:,:) = vn(:,:,:)                    ! (caution: the Asselin filter trends computation will be shifted by 1 timestep)
169      ENDIF
170
171      ! Time filter and swap of dynamics arrays
172      ! ---------------------------------------
173      IF( l_1st_euler ) THEN        !==  Euler at 1st time-step  ==!   (swap only)
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                          ! e3._n <-- e3._a
179            DO jk = 1, jpkm1
180               e3t_n(:,:,jk) = e3t_a(:,:,jk)
181               e3u_n(:,:,jk) = e3u_a(:,:,jk)
182               e3v_n(:,:,jk) = e3v_a(:,:,jk)
183            END DO
184         ENDIF
185         !
186      ELSE                          !==  Leap-Frog  ==!   (Asselin filter and swap)
187         !
188         !                                ! =============!
189         IF( ln_linssh ) THEN             ! Fixed volume !
190            !                             ! =============!
191            DO jk = 1, jpkm1                             
192               DO jj = 1, jpj
193                  DO ji = 1, jpi   
194                     zuf = un(ji,jj,jk) + rn_atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
195                     zvf = vn(ji,jj,jk) + rn_atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
196                     !
197                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
198                     vb(ji,jj,jk) = zvf
199                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
200                     vn(ji,jj,jk) = va(ji,jj,jk)
201                  END DO
202               END DO
203            END DO
204            !                             ! ================!
205         ELSE                             ! Variable volume !
206            !                             ! ================!
207            ! Before scale factor at t-points   (used as a now filtered scale factor until the swap)
208            DO jk = 1, jpkm1
209               e3t_b(:,:,jk) = e3t_n(:,:,jk) + rn_atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
210            END DO
211            ! Add volume filter correction: compatibility with tracer advection scheme
212            !    => time filter + conservation correction (only at the first level)
213            zcoef = rn_atfp * rn_Dt * r1_rho0
214
215            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
216
217            IF ( ln_rnf ) THEN
218               IF( ln_rnf_depth ) THEN
219                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
220                     DO jj = 1, jpj
221                        DO ji = 1, jpi
222                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
223                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
224                                  &              * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
225                           ENDIF
226                        END DO
227                     END DO
228                  END DO
229               ELSE
230                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
231               ENDIF
232            ENDIF
233
234            IF ( ln_isf ) THEN   ! if ice shelf melting
235               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
236                  DO jj = 1, jpj
237                     DO ji = 1, jpi
238                        IF( misfkt(ji,jj) <= jk .AND. jk <= misfkb(ji,jj) ) THEN
239                           e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
240                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
241                        ENDIF
242                     END DO
243                  END DO
244               END DO
245            ENDIF
246            !
247            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
248               ! Before filtered scale factor at (u/v)-points
249               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
250               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
251               DO jk = 1, jpkm1
252                  DO jj = 1, jpj
253                     DO ji = 1, jpi
254                        zuf = un(ji,jj,jk) + rn_atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
255                        zvf = vn(ji,jj,jk) + rn_atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
256                        !
257                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
258                        vb(ji,jj,jk) = zvf
259                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
260                        vn(ji,jj,jk) = va(ji,jj,jk)
261                     END DO
262                  END DO
263               END DO
264               !
265            ELSE                          ! Asselin filter applied on thickness weighted velocity
266               !
267               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
268               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
269               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
270               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
271               DO jk = 1, jpkm1
272                  DO jj = 1, jpj
273                     DO ji = 1, jpi                 
274                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
275                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
276                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
277                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
278                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
279                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
280                        !
281                        zuf = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
282                        zvf = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
283                        !
284                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
285                        vb(ji,jj,jk) = zvf
286                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
287                        vn(ji,jj,jk) = va(ji,jj,jk)
288                     END DO
289                  END DO
290               END DO
291               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
292               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
293               !
294               DEALLOCATE( ze3u_f , ze3v_f )
295            ENDIF
296            !
297         ENDIF
298         !
299         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
300            ! Revert "before" velocities to time split estimate
301            ! Doing it here also means that asselin filter contribution is removed 
302            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
303            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
304            DO jk = 2, jpkm1
305               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
306               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
307            END DO
308            DO jk = 1, jpkm1
309               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
310               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
311            END DO
312         ENDIF
313         !
314      ENDIF    ! end Leap-Frog time stepping
315      !
316      ! Set "now" and "before" barotropic velocities for next time step:
317      ! JC: Would be more clever to swap variables than to make a full vertical integration
318      !
319      !
320      IF(.NOT.ln_linssh ) THEN
321         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
322         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
323         DO jk = 2, jpkm1
324            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
325            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
326         END DO
327         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
328         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
329      ENDIF
330      !
331      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
332      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
333      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
334      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
335      DO jk = 2, jpkm1
336         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
337         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
338         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
339         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
340      END DO
341      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
342      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
343      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
344      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
345      !
346      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
347         CALL iom_put(  "ubar", un_b(:,:) )
348         CALL iom_put(  "vbar", vn_b(:,:) )
349      ENDIF
350      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
351         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * r1_Dt
352         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * r1_Dt
353         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
354      ENDIF
355      !
356      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
357         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
358      !
359      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
360      IF( l_trddyn     )   DEALLOCATE( zua, zva )
361      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
362      !
363   END SUBROUTINE dyn_nxt
364
365   !!=========================================================================
366END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.