source: branches/UKMO/GO6_dyn_vrt_diag/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 8300

Last change on this file since 8300 was 8300, checked in by glong, 4 years ago

Changed diagnostics to calculate the contributions by taking after-before before going into the specific numerical scheme for the model.

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