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/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 9150

Last change on this file since 9150 was 9150, checked in by deazer, 4 years ago

Sign fix and deal with bdy corner points.

  • Property svn:keywords set to Id
File size: 18.5 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 phycst         ! physical constants
32   USE dynadv         ! dynamics: vector invariant versus flux form
33   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
34   USE domvvl         ! variable volume
35   USE bdy_oce   , ONLY: ln_bdy
36   USE bdydta         ! ocean open boundary conditions
37   USE bdydyn         ! ocean open boundary conditions
38   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
39   USE trd_oce        ! trends: ocean variables
40   USE trddyn         ! trend manager: dynamics
41   USE trdken         ! trend manager: kinetic energy
42   !
43   USE in_out_manager ! I/O manager
44   USE iom            ! I/O manager library
45   USE lbclnk         ! lateral boundary condition (or mpp link)
46   USE lib_mpp        ! MPP library
47   USE wrk_nemo       ! Memory Allocation
48   USE prtctl         ! Print control
49   USE timing         ! Timing
50#if defined key_agrif
51   USE agrif_opa_interp
52#endif
53
54   IMPLICIT NONE
55   PRIVATE
56
57   PUBLIC    dyn_nxt   ! routine called by step.F90
58
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 (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) + 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, 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( l_trddyn           )   CALL wrk_alloc( jpi,jpj,jpk,   zua, zva)
108      !
109      IF( kt == nit000 ) THEN
110         IF(lwp) WRITE(numout,*)
111         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
112         IF(lwp) WRITE(numout,*) '~~~~~~~'
113      ENDIF
114
115      IF ( ln_dynspg_ts ) THEN
116         ! Ensure below that barotropic velocities match time splitting estimate
117         ! Compute actual transport and replace it with ts estimate at "after" time step
118         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
119         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
120         DO jk = 2, jpkm1
121            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
122            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
123         END DO
124         DO jk = 1, jpkm1
125            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
126            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
127         END DO
128         !
129         IF( .NOT.ln_bt_fw ) THEN
130            ! Remove advective velocity from "now velocities"
131            ! prior to asselin filtering     
132            ! In the forward case, this is done below after asselin filtering   
133            ! so that asselin contribution is removed at the same time
134            DO jk = 1, jpkm1
135               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
136               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
137            END DO 
138         ENDIF
139      ENDIF
140
141      ! Update after velocity on domain lateral boundaries
142      ! --------------------------------------------------     
143# if defined key_agrif
144      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
145# endif
146      !
147      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
148      CALL lbc_lnk( va, 'V', -1. ) 
149      !
150      !                                !* BDY open boundaries
151      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
152      IF( ln_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      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
157         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
158         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
159         !
160         !                                  ! Kinetic energy and Conversion
161         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
162         !
163         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
164            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
165            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
166            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
167            CALL iom_put( "vtrd_tot", zva )
168         ENDIF
169         !
170         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
171         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
172         !                                  !  computation of the asselin filter trends)
173      ENDIF
174
175      ! Time filter and swap of dynamics arrays
176      ! ------------------------------------------
177      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
178         DO jk = 1, jpkm1
179            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
180            vn(:,:,jk) = va(:,:,jk)
181         END DO
182         IF(.NOT.ln_linssh ) THEN
183            DO jk = 1, jpkm1
184               e3t_b(:,:,jk) = e3t_n(:,:,jk)
185               e3u_b(:,:,jk) = e3u_n(:,:,jk)
186               e3v_b(:,:,jk) = e3v_n(:,:,jk)
187            END DO
188         ENDIF
189      ELSE                                             !* Leap-Frog : Asselin filter and swap
190         !                                ! =============!
191         IF( ln_linssh ) THEN             ! Fixed volume !
192            !                             ! =============!
193            DO jk = 1, jpkm1                             
194               DO jj = 1, jpj
195                  DO ji = 1, jpi   
196                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
197                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
198                     !
199                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
200                     vb(ji,jj,jk) = zvf
201                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
202                     vn(ji,jj,jk) = va(ji,jj,jk)
203                  END DO
204               END DO
205            END DO
206            !                             ! ================!
207         ELSE                             ! Variable volume !
208            !                             ! ================!
209            ! Before scale factor at t-points
210            ! (used as a now filtered scale factor until the swap)
211            ! ----------------------------------------------------
212            IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN    ! No asselin filtering on thicknesses if forward time splitting
213               e3t_b(:,:,1:jpkm1) = e3t_n(:,:,1:jpkm1)
214            ELSE
215               DO jk = 1, jpkm1
216                  e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
217               END DO
218               ! Add volume filter correction: compatibility with tracer advection scheme
219               ! => time filter + conservation correction (only at the first level)
220               zcoef = atfp * rdt * r1_rau0
221               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting
222                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) &
223                                 &                      ) * tmask(:,:,1)
224                  IF( ln_rnf_depth ) THEN
225                     DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too, not sure for ice shelf case yet
226                        DO jj = 1, jpj
227                           DO ji = 1, jpi
228                              IF( mikt(ji,jj) <= jk .and. jk <=  nk_rnf(ji,jj)  ) THEN
229                                 e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef *  (-rnf_b(ji,jj) + rnf(ji,jj))*(e3t_n(ji,jj,jk)/h_rnf(ji,jj)  )*tmask(ji,jj,jk)
230                              ENDIF
231                           ENDDO
232                        ENDDO
233                     ENDDO
234                  ELSE
235                      e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  (-rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
236                  ENDIF
237               ELSE                     ! if ice shelf melting
238                  DO jj = 1, jpj
239                     DO ji = 1, jpi
240                        ikt = mikt(ji,jj)
241                        e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  &
242                           &                                           - rnf_b   (ji,jj) + rnf   (ji,jj)  &
243                           &                                           + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,ikt)
244                     END DO
245                  END DO
246               END IF
247            ENDIF
248            !
249            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
250               ! Before filtered scale factor at (u/v)-points
251               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
252               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
253               DO jk = 1, jpkm1
254                  DO jj = 1, jpj
255                     DO ji = 1, jpi
256                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
257                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
258                        !
259                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
260                        vb(ji,jj,jk) = zvf
261                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
262                        vn(ji,jj,jk) = va(ji,jj,jk)
263                     END DO
264                  END DO
265               END DO
266               !
267            ELSE                          ! Asselin filter applied on thickness weighted velocity
268               !
269               CALL wrk_alloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
270               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
271               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
272               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
273               DO jk = 1, jpkm1
274                  DO jj = 1, jpj
275                     DO ji = 1, jpi                 
276                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
277                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
278                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
279                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
280                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
281                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
282                        !
283                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
284                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
285                        !
286                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
287                        vb(ji,jj,jk) = zvf
288                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
289                        vn(ji,jj,jk) = va(ji,jj,jk)
290                     END DO
291                  END DO
292               END DO
293               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
294               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
295               !
296               CALL wrk_dealloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
297            ENDIF
298            !
299         ENDIF
300         !
301         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
302            ! Revert "before" velocities to time split estimate
303            ! Doing it here also means that asselin filter contribution is removed 
304            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
305            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
306            DO jk = 2, jpkm1
307               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
308               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
309            END DO
310            DO jk = 1, jpkm1
311               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
312               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
313            END DO
314         ENDIF
315         !
316      ENDIF ! neuler =/0
317      !
318      ! Set "now" and "before" barotropic velocities for next time step:
319      ! JC: Would be more clever to swap variables than to make a full vertical
320      ! integration
321      !
322      !
323      IF(.NOT.ln_linssh ) THEN
324         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
325         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
326         DO jk = 2, jpkm1
327            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
328            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
329         END DO
330         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
331         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
332      ENDIF
333      !
334      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
335      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
336      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
337      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
338      DO jk = 2, jpkm1
339         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
340         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
341         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
342         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
343      END DO
344      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
345      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
346      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
347      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
348      !
349      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
350         CALL iom_put(  "ubar", un_b(:,:) )
351         CALL iom_put(  "vbar", vn_b(:,:) )
352      ENDIF
353      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
354         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
355         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
356         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
357      ENDIF
358      !
359      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
360         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
361      !
362      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj,       zue, zve )
363      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk,   zua, zva )
364      !
365      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
366      !
367   END SUBROUTINE dyn_nxt
368
369   !!=========================================================================
370END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.