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

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 6006

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

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