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 NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert STO, TRD and USR modules and all knock on effects of these conversions. Note change to USR module may have implications for the TEST CASES (not tested yet). Standard SETTE tested only
  • 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.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 sbcisf         ! ice shelf
32   USE phycst         ! physical constants
33   USE dynadv         ! dynamics: vector invariant versus flux form
34   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
35   USE domvvl         ! variable volume
36   USE bdy_oce   , ONLY: ln_bdy
37   USE bdydta         ! ocean open boundary conditions
38   USE bdydyn         ! ocean open boundary conditions
39   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
40   USE trd_oce        ! trends: ocean variables
41   USE trddyn         ! trend manager: dynamics
42   USE trdken         ! trend manager: kinetic energy
43   !
44   USE in_out_manager ! I/O manager
45   USE iom            ! I/O manager library
46   USE lbclnk         ! lateral boundary condition (or mpp link)
47   USE lib_mpp        ! MPP library
48   USE prtctl         ! Print control
49   USE timing         ! Timing
50#if defined key_agrif
51   USE agrif_oce_interp
52#endif
53
54   IMPLICIT NONE
55   PRIVATE
56
57   PUBLIC    dyn_nxt   ! routine called by step.F90
58
59   !!----------------------------------------------------------------------
60   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
61   !! $Id$
62   !! Software governed by the CeCILL license (see ./LICENSE)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE dyn_nxt ( kt, Kmm )
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 (ln_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 non linear free surface,
88      !!             the time filter is applied on thickness weighted velocity.
89      !!             As a result, dyn_nxt MUST be called after tra_nxt.
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      INTEGER, INTENT( in ) ::   Kmm     ! time level index
96      !
97      INTEGER  ::   ji, jj, jk   ! dummy loop indices
98      INTEGER  ::   ikt          ! local integers
99      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
100      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
101      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
102      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
103      !!----------------------------------------------------------------------
104      !
105      IF( ln_timing    )   CALL timing_start('dyn_nxt')
106      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
107      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
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(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
119         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
120         DO jk = 2, jpkm1
121            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
122            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
123         END DO
124         DO jk = 1, jpkm1
125            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
126            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_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(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
136               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + 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_multi( 'dynnxt', ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
148      !
149      !                                !* BDY open boundaries
150      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
151      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
152
153!!$   Do we need a call to bdy_vol here??
154      !
155      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
156         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
157         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
158         !
159         !                                  ! Kinetic energy and Conversion
160         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt, Kmm )
161         !
162         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
163            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
164            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
165            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
166            CALL iom_put( "vtrd_tot", zva )
167         ENDIF
168         !
169         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
170         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
171         !                                  !  computation of the asselin filter trends)
172      ENDIF
173
174      ! Time filter and swap of dynamics arrays
175      ! ------------------------------------------
176      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
177         DO jk = 1, jpkm1
178            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
179            vn(:,:,jk) = va(:,:,jk)
180         END DO
181         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
182!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a 
183            DO jk = 1, jpkm1
184!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
185!               e3u_b(:,:,jk) = e3u_n(:,:,jk)
186!               e3v_b(:,:,jk) = e3v_n(:,:,jk)
187               !
188               e3t_n(:,:,jk) = e3t_a(:,:,jk)
189               e3u_n(:,:,jk) = e3u_a(:,:,jk)
190               e3v_n(:,:,jk) = e3v_a(:,:,jk)
191            END DO
192!!gm BUG end
193         ENDIF
194                                                            !
195         
196      ELSE                                             !* Leap-Frog : Asselin filter and swap
197         !                                ! =============!
198         IF( ln_linssh ) THEN             ! Fixed volume !
199            !                             ! =============!
200            DO jk = 1, jpkm1                             
201               DO jj = 1, jpj
202                  DO ji = 1, jpi   
203                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
204                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
205                     !
206                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
207                     vb(ji,jj,jk) = zvf
208                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
209                     vn(ji,jj,jk) = va(ji,jj,jk)
210                  END DO
211               END DO
212            END DO
213            !                             ! ================!
214         ELSE                             ! Variable volume !
215            !                             ! ================!
216            ! Before scale factor at t-points
217            ! (used as a now filtered scale factor until the swap)
218            ! ----------------------------------------------------
219            DO jk = 1, jpkm1
220               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
221            END DO
222            ! Add volume filter correction: compatibility with tracer advection scheme
223            ! => time filter + conservation correction (only at the first level)
224            zcoef = atfp * rdt * r1_rau0
225
226            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
227
228            IF ( ln_rnf ) THEN
229               IF( ln_rnf_depth ) THEN
230                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
231                     DO jj = 1, jpj
232                        DO ji = 1, jpi
233                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
234                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
235                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
236                           ENDIF
237                        ENDDO
238                     ENDDO
239                  ENDDO
240               ELSE
241                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
242               ENDIF
243            END IF
244
245            IF ( ln_isf ) THEN   ! if ice shelf melting
246               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
247                  DO jj = 1, jpj
248                     DO ji = 1, jpi
249                        IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN
250                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
251                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
252                        ELSEIF ( jk==misfkb(ji,jj) ) THEN
253                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
254                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
255                        ENDIF
256                     END DO
257                  END DO
258               END DO
259            END IF
260            !
261            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
262               ! Before filtered scale factor at (u/v)-points
263               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
264               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
265               DO jk = 1, jpkm1
266                  DO jj = 1, jpj
267                     DO ji = 1, jpi
268                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
269                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
270                        !
271                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
272                        vb(ji,jj,jk) = zvf
273                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
274                        vn(ji,jj,jk) = va(ji,jj,jk)
275                     END DO
276                  END DO
277               END DO
278               !
279            ELSE                          ! Asselin filter applied on thickness weighted velocity
280               !
281               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
282               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
283               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
284               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
285               DO jk = 1, jpkm1
286                  DO jj = 1, jpj
287                     DO ji = 1, jpi                 
288                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
289                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
290                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
291                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
292                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
293                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
294                        !
295                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
296                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
297                        !
298                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
299                        vb(ji,jj,jk) = zvf
300                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
301                        vn(ji,jj,jk) = va(ji,jj,jk)
302                     END DO
303                  END DO
304               END DO
305               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
306               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
307               !
308               DEALLOCATE( ze3u_f , ze3v_f )
309            ENDIF
310            !
311         ENDIF
312         !
313         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
314            ! Revert "before" velocities to time split estimate
315            ! Doing it here also means that asselin filter contribution is removed 
316            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
317            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
318            DO jk = 2, jpkm1
319               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
320               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
321            END DO
322            DO jk = 1, jpkm1
323               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
324               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
325            END DO
326         ENDIF
327         !
328      ENDIF ! neuler =/0
329      !
330      ! Set "now" and "before" barotropic velocities for next time step:
331      ! JC: Would be more clever to swap variables than to make a full vertical
332      ! integration
333      !
334      !
335      IF(.NOT.ln_linssh ) THEN
336         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
337         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
338         DO jk = 2, jpkm1
339            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
340            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
341         END DO
342         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
343         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
344      ENDIF
345      !
346      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
347      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
348      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
349      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
350      DO jk = 2, jpkm1
351         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
352         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
353         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
354         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
355      END DO
356      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
357      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
358      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
359      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
360      !
361      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
362         CALL iom_put(  "ubar", un_b(:,:) )
363         CALL iom_put(  "vbar", vn_b(:,:) )
364      ENDIF
365      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
366         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
367         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
368         CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm )
369      ENDIF
370      !
371      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
372         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
373      !
374      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
375      IF( l_trddyn     )   DEALLOCATE( zua, zva )
376      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
377      !
378   END SUBROUTINE dyn_nxt
379
380   !!=========================================================================
381END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.