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

source: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 10684

Last change on this file since 10684 was 10684, checked in by jcastill, 6 years ago

Remove svn keywords
File size: 17.7 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 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 (ln_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      !                                !* BDY open boundaries
150      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
151      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
152
153!!$   Do we need a call to bdy_vol here??
154      !
155      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
156         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
157         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
158         !
159         !                                  ! Kinetic energy and Conversion
160         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
161         !
162         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
163            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
164            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
165            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
166            CALL iom_put( "vtrd_tot", zva )
167         ENDIF
168         !
169         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
170         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
171         !                                  !  computation of the asselin filter trends)
172      ENDIF
173
174      ! Time filter and swap of dynamics arrays
175      ! ------------------------------------------
176      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
177         DO jk = 1, jpkm1
178            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
179            vn(:,:,jk) = va(:,:,jk)
180         END DO
181         IF(.NOT.ln_linssh ) THEN
182            DO jk = 1, jpkm1
183               e3t_b(:,:,jk) = e3t_n(:,:,jk)
184               e3u_b(:,:,jk) = e3u_n(:,:,jk)
185               e3v_b(:,:,jk) = e3v_n(:,:,jk)
186            END DO
187         ENDIF
188      ELSE                                             !* Leap-Frog : Asselin filter and swap
189         !                                ! =============!
190         IF( ln_linssh ) THEN             ! Fixed volume !
191            !                             ! =============!
192            DO jk = 1, jpkm1                             
193               DO jj = 1, jpj
194                  DO ji = 1, jpi   
195                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
196                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
197                     !
198                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
199                     vb(ji,jj,jk) = zvf
200                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
201                     vn(ji,jj,jk) = va(ji,jj,jk)
202                  END DO
203               END DO
204            END DO
205            !                             ! ================!
206         ELSE                             ! Variable volume !
207            !                             ! ================!
208            ! Before scale factor at t-points
209            ! (used as a now filtered scale factor until the swap)
210            ! ----------------------------------------------------
211            IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN    ! No asselin filtering on thicknesses if forward time splitting
212               e3t_b(:,:,1:jpkm1) = e3t_n(:,:,1:jpkm1)
213            ELSE
214               DO jk = 1, jpkm1
215                  e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
216               END DO
217               ! Add volume filter correction: compatibility with tracer advection scheme
218               ! => time filter + conservation correction (only at the first level)
219               zcoef = atfp * rdt * r1_rau0
220               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting
221                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) &
222                                 &                      - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
223               ELSE                     ! if ice shelf melting
224                  DO jj = 1, jpj
225                     DO ji = 1, jpi
226                        ikt = mikt(ji,jj)
227                        e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  &
228                           &                                           - rnf_b   (ji,jj) + rnf   (ji,jj)  &
229                           &                                           + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,ikt)
230                     END DO
231                  END DO
232               END IF
233            ENDIF
234            !
235            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
236               ! Before filtered scale factor at (u/v)-points
237               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
238               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
239               DO jk = 1, jpkm1
240                  DO jj = 1, jpj
241                     DO ji = 1, jpi
242                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
243                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
244                        !
245                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
246                        vb(ji,jj,jk) = zvf
247                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
248                        vn(ji,jj,jk) = va(ji,jj,jk)
249                     END DO
250                  END DO
251               END DO
252               !
253            ELSE                          ! Asselin filter applied on thickness weighted velocity
254               !
255               CALL wrk_alloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
256               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
257               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
258               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
259               DO jk = 1, jpkm1
260                  DO jj = 1, jpj
261                     DO ji = 1, jpi                 
262                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
263                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
264                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
265                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
266                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
267                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
268                        !
269                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
270                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
271                        !
272                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
273                        vb(ji,jj,jk) = zvf
274                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
275                        vn(ji,jj,jk) = va(ji,jj,jk)
276                     END DO
277                  END DO
278               END DO
279               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
280               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
281               !
282               CALL wrk_dealloc( jpi,jpj,jpk,   ze3u_f, ze3v_f )
283            ENDIF
284            !
285         ENDIF
286         !
287         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
288            ! Revert "before" velocities to time split estimate
289            ! Doing it here also means that asselin filter contribution is removed 
290            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
291            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
292            DO jk = 2, jpkm1
293               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
294               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
295            END DO
296            DO jk = 1, jpkm1
297               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
298               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
299            END DO
300         ENDIF
301         !
302      ENDIF ! neuler =/0
303      !
304      ! Set "now" and "before" barotropic velocities for next time step:
305      ! JC: Would be more clever to swap variables than to make a full vertical
306      ! integration
307      !
308      !
309      IF(.NOT.ln_linssh ) THEN
310         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
311         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
312         DO jk = 2, jpkm1
313            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
314            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
315         END DO
316         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
317         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
318      ENDIF
319      !
320      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
321      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
322      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
323      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
324      DO jk = 2, jpkm1
325         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
326         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
327         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
328         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
329      END DO
330      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
331      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
332      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
333      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
334      !
335      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
336         CALL iom_put(  "ubar", un_b(:,:) )
337         CALL iom_put(  "vbar", vn_b(:,:) )
338      ENDIF
339      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
340         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
341         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
342         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
343      ENDIF
344      !
345      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
346         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
347      !
348      IF( ln_dynspg_ts )   CALL wrk_dealloc( jpi,jpj,       zue, zve )
349      IF( l_trddyn     )   CALL wrk_dealloc( jpi,jpj,jpk,   zua, zva )
350      !
351      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
352      !
353   END SUBROUTINE dyn_nxt
354
355   !!=========================================================================
356END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.