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

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 5945

Last change on this file since 5945 was 5945, checked in by mathiot, 8 years ago

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