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

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 9294

Last change on this file since 9294 was 5628, checked in by mathiot, 9 years ago

correction of unconservation with ice shelf following corrections did by Jerome on the runoff

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