source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 9119

Last change on this file since 9119 was 9119, checked in by nicolasmartin, 3 years ago

Fix longer lines so should be harmless (passed SETTE compilations)

  • Property svn:keywords set to Id
File size: 19.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.6  !  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 sbcrnf         ! river runoffs
31   USE phycst         ! physical constants
32   USE dynadv         ! dynamics: vector invariant versus flux form
33   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
34   USE domvvl         ! variable volume
35   USE bdy_oce   , ONLY: ln_bdy
36   USE bdydta         ! ocean open boundary conditions
37   USE bdydyn         ! ocean open boundary conditions
38   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
39   USE trd_oce        ! trends: ocean variables
40   USE trddyn         ! trend manager: dynamics
41   USE trdken         ! trend manager: kinetic energy
42   !
43   USE in_out_manager ! I/O manager
44   USE iom            ! I/O manager library
45   USE lbclnk         ! lateral boundary condition (or mpp link)
46   USE lib_mpp        ! MPP library
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   !!----------------------------------------------------------------------
59   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
60   !! $Id$
61   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   SUBROUTINE dyn_nxt ( kt )
66      !!----------------------------------------------------------------------
67      !!                  ***  ROUTINE dyn_nxt  ***
68      !!                   
69      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
70      !!             condition on the after velocity, achieve the time stepping
71      !!             by applying the Asselin filter on now fields and swapping
72      !!             the fields.
73      !!
74      !! ** Method  : * Ensure after velocities transport matches time splitting
75      !!             estimate (ln_dynspg_ts=T)
76      !!
77      !!              * Apply lateral boundary conditions on after velocity
78      !!             at the local domain boundaries through lbc_lnk call,
79      !!             at the one-way open boundaries (ln_bdy=T),
80      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
81      !!
82      !!              * Apply the time filter applied and swap of the dynamics
83      !!             arrays to start the next time step:
84      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
85      !!                (un,vn) = (ua,va).
86      !!             Note that with flux form advection and non linear free surface,
87      !!             the time filter is applied on thickness weighted velocity.
88      !!             As a result, dyn_nxt MUST be called after tra_nxt.
89      !!
90      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
91      !!               un,vn   now horizontal velocity of next time-step
92      !!----------------------------------------------------------------------
93      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
94      !
95      INTEGER  ::   ji, jj, jk   ! dummy loop indices
96      INTEGER  ::   ikt          ! local integers
97      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
98      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
99      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
100      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
101      !!----------------------------------------------------------------------
102      !
103      IF( ln_timing    )   CALL timing_start('dyn_nxt')
104      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
105      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
106      !
107      IF( kt == nit000 ) THEN
108         IF(lwp) WRITE(numout,*)
109         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
110         IF(lwp) WRITE(numout,*) '~~~~~~~'
111      ENDIF
112
113      IF ( ln_dynspg_ts ) THEN
114         ! Ensure below that barotropic velocities match time splitting estimate
115         ! Compute actual transport and replace it with ts estimate at "after" time step
116         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
117         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
118         DO jk = 2, jpkm1
119            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
120            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
121         END DO
122         DO jk = 1, jpkm1
123            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
124            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
125         END DO
126         !
127         IF( .NOT.ln_bt_fw ) THEN
128            ! Remove advective velocity from "now velocities"
129            ! prior to asselin filtering     
130            ! In the forward case, this is done below after asselin filtering   
131            ! so that asselin contribution is removed at the same time
132            DO jk = 1, jpkm1
133               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
134               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
135            END DO 
136         ENDIF
137      ENDIF
138
139      ! Update after velocity on domain lateral boundaries
140      ! --------------------------------------------------     
141# if defined key_agrif
142      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
143# endif
144      !
145      CALL lbc_lnk_multi( ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
146      !
147      !                                !* BDY open boundaries
148      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
149      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
150
151!!$   Do we need a call to bdy_vol here??
152      !
153      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
154         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
155         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
156         !
157         !                                  ! Kinetic energy and Conversion
158         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
159         !
160         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
161            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
162            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
163            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
164            CALL iom_put( "vtrd_tot", zva )
165         ENDIF
166         !
167         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
168         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
169         !                                  !  computation of the asselin filter trends)
170      ENDIF
171
172      ! Time filter and swap of dynamics arrays
173      ! ------------------------------------------
174      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
175         DO jk = 1, jpkm1
176            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
177            vn(:,:,jk) = va(:,:,jk)
178         END DO
179         IF(.NOT.ln_linssh ) THEN
180            DO jk = 1, jpkm1
181               e3t_b(:,:,jk) = e3t_n(:,:,jk)
182               e3u_b(:,:,jk) = e3u_n(:,:,jk)
183               e3v_b(:,:,jk) = e3v_n(:,:,jk)
184            END DO
185         ENDIF
186      ELSE                                             !* Leap-Frog : Asselin filter and swap
187         !                                ! =============!
188         IF( ln_linssh ) THEN             ! Fixed volume !
189            !                             ! =============!
190            DO jk = 1, jpkm1                             
191               DO jj = 1, jpj
192                  DO ji = 1, jpi   
193                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
194                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
195                     !
196                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
197                     vb(ji,jj,jk) = zvf
198                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
199                     vn(ji,jj,jk) = va(ji,jj,jk)
200                  END DO
201               END DO
202            END DO
203            !                             ! ================!
204         ELSE                             ! Variable volume !
205            !                             ! ================!
206            ! Before scale factor at t-points
207            ! (used as a now filtered scale factor until the swap)
208            ! ----------------------------------------------------
209            DO jk = 1, jpkm1
210               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
211            END DO
212            ! Add volume filter correction: compatibility with tracer advection scheme
213            ! => time filter + conservation correction (only at the first level)
214            zcoef = atfp * rdt * r1_rau0
215            IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting
216               e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) &
217                              &                      - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
218            ELSE                     ! if ice shelf melting
219               DO jj = 1, jpj
220                  DO ji = 1, jpi
221                     ikt = mikt(ji,jj)
222                     e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  &
223                        &                                           - rnf_b   (ji,jj) + rnf   (ji,jj)  &
224                        &                                           + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,ikt)
225                  END DO
226               END DO
227               ! Add volume filter correction: compatibility with tracer advection scheme
228               ! => time filter + conservation correction (only at the first level)
229               zcoef = atfp * rdt * r1_rau0
230               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting
231                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) &
232                                 &                      ) * tmask(:,:,1)
233                  IF( ln_rnf_depth ) THEN
234                     DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too, not sure for ice shelf case yet
235                        DO jj = 1, jpj
236                           DO ji = 1, jpi
237                              IF( mikt(ji,jj) <= jk .and. jk <=  nk_rnf(ji,jj)  ) THEN
238                                 e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
239                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
240                              ENDIF
241                           ENDDO
242                        ENDDO
243                     ENDDO
244                  ELSE
245                      e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
246                  ENDIF
247               ELSE                     ! if ice shelf melting
248                  DO jj = 1, jpj
249                     DO ji = 1, jpi
250                        ikt = mikt(ji,jj)
251                        e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  &
252                           &                                           - rnf_b   (ji,jj) + rnf   (ji,jj)  &
253                           &                                           + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,ikt)
254                     END DO
255                  END DO
256               END IF
257            END IF
258            !
259            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
260               ! Before filtered scale factor at (u/v)-points
261               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
262               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
263               DO jk = 1, jpkm1
264                  DO jj = 1, jpj
265                     DO ji = 1, jpi
266                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
267                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
268                        !
269                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
270                        vb(ji,jj,jk) = zvf
271                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
272                        vn(ji,jj,jk) = va(ji,jj,jk)
273                     END DO
274                  END DO
275               END DO
276               !
277            ELSE                          ! Asselin filter applied on thickness weighted velocity
278               !
279               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
280               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
281               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
282               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
283               DO jk = 1, jpkm1
284                  DO jj = 1, jpj
285                     DO ji = 1, jpi                 
286                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
287                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
288                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
289                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
290                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
291                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
292                        !
293                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
294                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
295                        !
296                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
297                        vb(ji,jj,jk) = zvf
298                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
299                        vn(ji,jj,jk) = va(ji,jj,jk)
300                     END DO
301                  END DO
302               END DO
303               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
304               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
305               !
306               DEALLOCATE( ze3u_f , ze3v_f )
307            ENDIF
308            !
309         ENDIF
310         !
311         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
312            ! Revert "before" velocities to time split estimate
313            ! Doing it here also means that asselin filter contribution is removed 
314            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
315            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
316            DO jk = 2, jpkm1
317               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
318               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
319            END DO
320            DO jk = 1, jpkm1
321               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
322               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
323            END DO
324         ENDIF
325         !
326      ENDIF ! neuler =/0
327      !
328      ! Set "now" and "before" barotropic velocities for next time step:
329      ! JC: Would be more clever to swap variables than to make a full vertical
330      ! integration
331      !
332      !
333      IF(.NOT.ln_linssh ) THEN
334         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
335         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
336         DO jk = 2, jpkm1
337            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
338            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
339         END DO
340         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
341         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
342      ENDIF
343      !
344      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
345      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
346      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
347      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
348      DO jk = 2, jpkm1
349         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
350         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
351         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
352         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
353      END DO
354      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
355      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
356      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
357      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
358      !
359      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
360         CALL iom_put(  "ubar", un_b(:,:) )
361         CALL iom_put(  "vbar", vn_b(:,:) )
362      ENDIF
363      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
364         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
365         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
366         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
367      ENDIF
368      !
369      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
370         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
371      !
372      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
373      IF( l_trddyn     )   DEALLOCATE( zua, zva )
374      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
375      !
376   END SUBROUTINE dyn_nxt
377
378   !!=========================================================================
379END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.