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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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