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 on Ticket #1584 – Attachment – NEMO

Ticket #1584: dynnxt.F90

File dynnxt.F90, 19.0 KB (added by Robin_Waldman, 5 years ago)
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.6  !  2019-09  (R. Waldman) remove the barotropic trend update from time splitting  (added to dynzdf.F90)
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 dynspg_oce      ! type of surface pressure gradient
32   USE dynadv          ! dynamics: vector invariant versus flux form
33   USE domvvl          ! variable volume
34   USE bdy_oce         ! ocean open boundary conditions
35   USE bdydta          ! ocean open boundary conditions
36   USE bdydyn          ! ocean open boundary conditions
37   USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
38   USE trd_oce         ! trends: ocean variables
39   USE trddyn          ! trend manager: dynamics
40   USE trdken          ! trend manager: kinetic energy
41   !
42   USE in_out_manager  ! I/O manager
43   USE iom             ! I/O manager library
44   USE lbclnk          ! lateral boundary condition (or mpp link)
45   USE lib_mpp         ! MPP library
46   USE wrk_nemo        ! Memory Allocation
47   USE prtctl          ! Print control
48   USE timing          ! Timing
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      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
110      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
111      !!----------------------------------------------------------------------
112      !
113      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
114      !
115      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
116      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
117      !
118      IF( kt == nit000 ) THEN
119         IF(lwp) WRITE(numout,*)
120         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
121         IF(lwp) WRITE(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      IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN
161         ! Remove advective velocity from "now velocities"
162         ! prior to asselin filtering     
163         ! In the forward case, this is done below after asselin filtering   
164         ! so that asselin contribution is removed at the same time
165         DO jk = 1, jpkm1
166            un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
167            vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
168         END DO 
169      ENDIF
170# endif
171
172      ! Update after velocity on domain lateral boundaries
173      ! --------------------------------------------------     
174      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
175      CALL lbc_lnk( va, 'V', -1. ) 
176      !
177# if defined key_bdy
178      !                                !* BDY open boundaries
179      IF( lk_bdy .AND. lk_dynspg_exp ) CALL bdy_dyn( kt )
180      IF( lk_bdy .AND. lk_dynspg_ts  ) CALL bdy_dyn( kt, dyn3d_only=.true. )
181
182!!$   Do we need a call to bdy_vol here??
183      !
184# endif
185      !
186# if defined key_agrif
187      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
188# endif
189#endif
190
191      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
192         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
193         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
194         !
195         !                                  ! Kinetic energy and Conversion
196         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
197         !
198         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
199            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
200            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
201            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
202            CALL iom_put( "vtrd_tot", zva )
203         ENDIF
204         !
205         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
206         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
207         !                                  !  computation of the asselin filter trends)
208      ENDIF
209
210      ! Time filter and swap of dynamics arrays
211      ! ------------------------------------------
212      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
213         DO jk = 1, jpkm1
214            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
215            vn(:,:,jk) = va(:,:,jk)
216         END DO
217         IF (lk_vvl) THEN
218            DO jk = 1, jpkm1
219               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
220               fse3u_b(:,:,jk) = fse3u_n(:,:,jk)
221               fse3v_b(:,:,jk) = fse3v_n(:,:,jk)
222            ENDDO
223         ENDIF
224      ELSE                                             !* Leap-Frog : Asselin filter and swap
225         !                                ! =============!
226         IF( .NOT. lk_vvl ) THEN          ! Fixed volume !
227            !                             ! =============!
228            DO jk = 1, jpkm1                             
229               DO jj = 1, jpj
230                  DO ji = 1, jpi   
231                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
232                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
233                     !
234                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
235                     vb(ji,jj,jk) = zvf
236                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
237                     vn(ji,jj,jk) = va(ji,jj,jk)
238                  END DO
239               END DO
240            END DO
241            !                             ! ================!
242         ELSE                             ! Variable volume !
243            !                             ! ================!
244            ! Before scale factor at t-points
245            ! (used as a now filtered scale factor until the swap)
246            ! ----------------------------------------------------
247            IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
248               ! No asselin filtering on thicknesses if forward time splitting
249                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
250            ELSE
251               fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) )
252               ! Add volume filter correction: compatibility with tracer advection scheme
253               ! => time filter + conservation correction (only at the first level)
254               IF ( nn_isf == 0) THEN   ! if no ice shelf melting
255                  fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) &
256                                 &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
257               ELSE                     ! if ice shelf melting
258                  DO jj = 1,jpj
259                     DO ji = 1,jpi
260                        jk = mikt(ji,jj)
261                        fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       &
262                                          &                          * ( (emp_b(ji,jj)    - emp(ji,jj)   ) &
263                                          &                            - (rnf_b(ji,jj)    - rnf(ji,jj)   ) &
264                                          &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk)
265                     END DO
266                  END DO
267               END IF
268            ENDIF
269            !
270            IF( ln_dynadv_vec ) THEN
271               ! Before scale factor at (u/v)-points
272               ! -----------------------------------
273               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
274               CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
275               ! Leap-Frog - Asselin filter and swap: applied on velocity
276               ! -----------------------------------
277               DO jk = 1, jpkm1
278                  DO jj = 1, jpj
279                     DO ji = 1, jpi
280                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
281                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
282                        !
283                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
284                        vb(ji,jj,jk) = zvf
285                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
286                        vn(ji,jj,jk) = va(ji,jj,jk)
287                     END DO
288                  END DO
289               END DO
290               !
291            ELSE
292               ! Temporary filtered scale factor at (u/v)-points (will become before scale factor)
293               !------------------------------------------------
294               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' )
295               CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' )
296               ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity
297               ! -----------------------------------             ===========================
298               DO jk = 1, jpkm1
299                  DO jj = 1, jpj
300                     DO ji = 1, jpi                 
301                        zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk)
302                        zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk)
303                        zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk)
304                        zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk)
305                        zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk)
306                        zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk)
307                        !
308                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
309                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
310                        !
311                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
312                        vb(ji,jj,jk) = zvf
313                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
314                        vn(ji,jj,jk) = va(ji,jj,jk)
315                     END DO
316                  END DO
317               END DO
318               fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)      ! e3u_b <-- filtered scale factor
319               fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
320            ENDIF
321            !
322         ENDIF
323         !
324         IF (lk_dynspg_ts.AND.ln_bt_fw) THEN
325            ! Revert "before" velocities to time split estimate
326            ! Doing it here also means that asselin filter contribution is removed 
327            zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
328            zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
329            DO jk = 2, jpkm1
330               zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
331               zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
332            END DO
333            DO jk = 1, jpkm1
334               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
335               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
336            END DO
337         ENDIF
338         !
339      ENDIF ! neuler =/0
340      !
341      ! Set "now" and "before" barotropic velocities for next time step:
342      ! JC: Would be more clever to swap variables than to make a full vertical
343      ! integration
344      !
345      !
346      IF (lk_vvl) THEN
347         hu_b(:,:) = 0.
348         hv_b(:,:) = 0.
349         DO jk = 1, jpkm1
350            hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
351            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
352         END DO
353         hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
354         hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
355      ENDIF
356      !
357      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
358      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
359      !
360      DO jk = 1, jpkm1
361         DO jj = 1, jpj
362            DO ji = 1, jpi
363               un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
364               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
365               !
366               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
367               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
368            END DO
369         END DO
370      END DO
371      !
372      !
373      un_b(:,:) = un_b(:,:) * hur_a(:,:)
374      vn_b(:,:) = vn_b(:,:) * hvr_a(:,:)
375      ub_b(:,:) = ub_b(:,:) * hur_b(:,:)
376      vb_b(:,:) = vb_b(:,:) * hvr_b(:,:)
377      !
378      !
379
380      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
381         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
382         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
383         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
384      ENDIF
385      !
386      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
387         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
388      !
389      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
390      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
391      !
392      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
393      !
394   END SUBROUTINE dyn_nxt
395
396   !!=========================================================================
397END MODULE dynnxt