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

source: branches/2015/dev_r5847_MERCATOR9_solveur_simplification/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 5868

Last change on this file since 5868 was 5868, checked in by jchanut, 8 years ago

Free surface simplification #1620. Step 1: suppress filtered free surface

  • Property svn:keywords set to Id
File size: 19.0 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.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
21   !!-------------------------------------------------------------------------
22 
23   !!-------------------------------------------------------------------------
24   !!   dyn_nxt      : obtain the next (after) horizontal velocity
25   !!-------------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE sbc_oce         ! Surface boundary condition: ocean fields
29   USE phycst          ! physical constants
30   USE dynspg_oce      ! type of surface pressure gradient
31   USE dynadv          ! dynamics: vector invariant versus flux form
32   USE domvvl          ! variable volume
33   USE bdy_oce         ! ocean open boundary conditions
34   USE bdydta          ! ocean open boundary conditions
35   USE bdydyn          ! ocean open boundary conditions
36   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
37   USE trd_oce         ! trends: ocean variables
38   USE trddyn          ! trend manager: dynamics
39   USE trdken          ! trend manager: kinetic energy
40   !
41   USE in_out_manager  ! I/O manager
42   USE iom             ! I/O manager library
43   USE lbclnk          ! lateral boundary condition (or mpp link)
44   USE lib_mpp         ! MPP library
45   USE wrk_nemo        ! Memory Allocation
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   !! * Substitutions
58#  include "domzgr_substitute.h90"
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 :   Compute the after horizontal velocity. Apply the boundary
71      !!             condition on the after velocity, achieved the time stepping
72      !!             by applying the Asselin filter on now fields and swapping
73      !!             the fields.
74      !!
75      !! ** Method  : * After velocity is compute using a leap-frog scheme:
76      !!                       (ua,va) = (ub,vb) + 2 rdt (ua,va)
77      !!             Note that with flux form advection and variable volume layer
78      !!             (lk_vvl=T), the leap-frog is applied on thickness weighted
79      !!             velocity.
80      !!
81      !!              * Apply lateral boundary conditions on after velocity
82      !!             at the local domain boundaries through lbc_lnk call,
83      !!             at the one-way open boundaries (lk_bdy=T),
84      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
85      !!
86      !!              * Apply the time filter applied and swap of the dynamics
87      !!             arrays to start the next time step:
88      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
89      !!                (un,vn) = (ua,va).
90      !!             Note that with flux form advection and variable volume layer
91      !!             (lk_vvl=T), the time filter is applied on thickness weighted
92      !!             velocity.
93      !!
94      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
95      !!               un,vn   now horizontal velocity of next time-step
96      !!----------------------------------------------------------------------
97      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
98      !
99      INTEGER  ::   ji, jj, jk   ! dummy loop indices
100      INTEGER  ::   iku, ikv     ! local integers
101      REAL(wp) ::   z2dt         ! temporary scalar
102      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars
103      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
104      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
105      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
106      !!----------------------------------------------------------------------
107      !
108      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
109      !
110      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
111      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
112      !
113      IF( kt == nit000 ) THEN
114         IF(lwp) WRITE(numout,*)
115         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
116         IF(lwp) WRITE(numout,*) '~~~~~~~'
117      ENDIF
118
119# if defined key_dynspg_exp
120      ! Next velocity :   Leap-frog time stepping
121      ! -------------
122      z2dt = 2. * rdt                                 ! Euler or leap-frog time step
123      IF( neuler == 0 .AND. kt == nit000 )  z2dt = rdt
124      !
125      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
126         DO jk = 1, jpkm1
127            ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
128            va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
129         END DO
130      ELSE                                            ! applied on thickness weighted velocity
131         DO jk = 1, jpkm1
132            ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)      &
133               &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk)  )   &
134               &         / fse3u_a(:,:,jk) * umask(:,:,jk)
135            va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)      &
136               &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk)  )   &
137               &         / fse3v_a(:,:,jk) * vmask(:,:,jk)
138         END DO
139      ENDIF
140# endif
141
142# if defined key_dynspg_ts
143!!gm IF ( lk_dynspg_ts ) THEN ....
144      ! Ensure below that barotropic velocities match time splitting estimate
145      ! Compute actual transport and replace it with ts estimate at "after" time step
146      zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
147      zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
148      DO jk = 2, jpkm1
149         zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
150         zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
151      END DO
152      DO jk = 1, jpkm1
153         ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
154         va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
155      END DO
156
157      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
158         ! Remove advective velocity from "now velocities"
159         ! prior to asselin filtering     
160         ! In the forward case, this is done below after asselin filtering   
161         ! so that asselin contribution is removed at the same time
162         DO jk = 1, jpkm1
163            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
164            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
165         END DO 
166      ENDIF
167!!gm ENDIF
168# endif
169
170      ! Update after velocity on domain lateral boundaries
171      ! --------------------------------------------------     
172      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
173      CALL lbc_lnk( va, 'V', -1. ) 
174      !
175# if defined key_bdy
176      !                                !* BDY open boundaries
177      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
178      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
179
180!!$   Do we need a call to bdy_vol here??
181      !
182# endif
183      !
184# if defined key_agrif
185      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
186# endif
187
188      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
189         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
190         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
191         !
192         !                                  ! Kinetic energy and Conversion
193         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
194         !
195         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
196            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
197            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
198            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
199            CALL iom_put( "vtrd_tot", zva )
200         ENDIF
201         !
202         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
203         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
204         !                                  !  computation of the asselin filter trends)
205      ENDIF
206
207      ! Time filter and swap of dynamics arrays
208      ! ------------------------------------------
209      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
210         DO jk = 1, jpkm1
211            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
212            vn(:,:,jk) = va(:,:,jk)
213         END DO
214         IF (lk_vvl) THEN
215            DO jk = 1, jpkm1
216               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
217               fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
218               fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
219            ENDDO
220         ENDIF
221      ELSE                                             !* Leap-Frog : Asselin filter and swap
222         !                                ! =============!
223         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
224            !                             ! =============!
225            DO jk = 1, jpkm1                             
226               DO jj = 1, jpj
227                  DO ji = 1, jpi   
228                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
229                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
230                     !
231                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
232                     vb(ji,jj,jk) = zvf
233                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
234                     vn(ji,jj,jk) = va(ji,jj,jk)
235                  END DO
236               END DO
237            END DO
238            !                             ! ================!
239         ELSE                             ! Variable volume !
240            !                             ! ================!
241            ! Before scale factor at t-points
242            ! (used as a now filtered scale factor until the swap)
243            ! ----------------------------------------------------
244            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
245               ! No asselin filtering on thicknesses if forward time splitting
246                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
247            ELSE
248               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
249               ! Add volume filter correction: compatibility with tracer advection scheme
250               ! => time filter + conservation correction (only at the first level)
251               IF ( nn_isf == 0) THEN   ! if no ice shelf melting
252                  fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
253                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
254               ELSE                     ! if ice shelf melting
255                  DO jj = 1,jpj
256                     DO ji = 1,jpi
257                        jk = mikt(ji,jj)
258                        fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       &
259                                          &                          * ( (emp_b(ji,jj)    - emp(ji,jj)   ) &
260                                          &                            - (rnf_b(ji,jj)    - rnf(ji,jj)   ) &
261                                          &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk)
262                     END DO
263                  END DO
264               END IF
265            ENDIF
266            !
267            IF( ln_dynadv_vec ) THEN
268               ! Before scale factor at (u/v)-points
269               ! -----------------------------------
270               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
271               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
272               ! Leap-Frog - Asselin filter and swap: applied on velocity
273               ! -----------------------------------
274               DO jk = 1, jpkm1
275                  DO jj = 1, jpj
276                     DO ji = 1, jpi
277                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
278                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
279                        !
280                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
281                        vb(ji,jj,jk) = zvf
282                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
283                        vn(ji,jj,jk) = va(ji,jj,jk)
284                     END DO
285                  END DO
286               END DO
287               !
288            ELSE
289               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
290               !------------------------------------------------
291               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
292               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
293               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
294               ! -----------------------------------             ===========================
295               DO jk = 1, jpkm1
296                  DO jj = 1, jpj
297                     DO ji = 1, jpi                 
298                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
299                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
300                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
301                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
302                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
303                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
304                        !
305                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
306                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
307                        !
308                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
309                        vb(ji,jj,jk) = zvf
310                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
311                        vn(ji,jj,jk) = va(ji,jj,jk)
312                     END DO
313                  END DO
314               END DO
315               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
316               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
317            ENDIF
318            !
319         ENDIF
320         !
321         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
322            ! Revert "before" velocities to time split estimate
323            ! Doing it here also means that asselin filter contribution is removed 
324            zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
325            zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
326            DO jk = 2, jpkm1
327               zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
328               zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
329            END DO
330            DO jk = 1, jpkm1
331               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
332               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
333            END DO
334         ENDIF
335         !
336      ENDIF ! neuler =/0
337      !
338      ! Set "now" and "before" barotropic velocities for next time step:
339      ! JC: Would be more clever to swap variables than to make a full vertical
340      ! integration
341      !
342      !
343      IF (lk_vvl) THEN
344         hu_b(:,:) = 0.
345         hv_b(:,:) = 0.
346         DO jk = 1, jpkm1
347            hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
348            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
349         END DO
350         hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
351         hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
352      ENDIF
353      !
354      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
355      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
356      !
357      DO jk = 1, jpkm1
358         DO jj = 1, jpj
359            DO ji = 1, jpi
360               un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
361               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
362               !
363               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
364               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
365            END DO
366         END DO
367      END DO
368      !
369      !
370      un_b(:,:) = un_b(:,:) * hur_a(:,:)
371      vn_b(:,:) = vn_b(:,:) * hvr_a(:,:)
372      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
373      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
374      !
375      !
376
377      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
378         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
379         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
380         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
381      ENDIF
382      !
383      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
384         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
385      !
386      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
387      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
388      !
389      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
390      !
391   END SUBROUTINE dyn_nxt
392
393   !!=========================================================================
394END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.