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 @ 10957

Last change on this file since 10957 was 10957, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : BDY module (including some removal of redundant code). Tested in AMM12 and ORCA1 (not full SETTE test at this stage).
  • Property svn:keywords set to Id
File size: 19.3 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, Kbb, Kmm, Kaa )
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 ) ::   Kbb, Kmm, Kaa  ! time level indices
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      !! IMMERSE development : Following the general pattern for the new code we want to pass in the
151      !!                       velocities to bdy_dyn as arguments so here we use "uu" and "vv" even
152      !!                       though we haven't converted the velocity names in the rest of dynnxt.F90
153      !!                       because it will be completely rewritten. DS.
154      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, uu, vv, Kaa )
155      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, uu, vv, Kaa, dyn3d_only=.true. )
156
157!!$   Do we need a call to bdy_vol here??
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, Kmm )
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( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
186!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a 
187            DO jk = 1, jpkm1
188!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
189!               e3u_b(:,:,jk) = e3u_n(:,:,jk)
190!               e3v_b(:,:,jk) = e3v_n(:,:,jk)
191               !
192               e3t_n(:,:,jk) = e3t_a(:,:,jk)
193               e3u_n(:,:,jk) = e3u_a(:,:,jk)
194               e3v_n(:,:,jk) = e3v_a(:,:,jk)
195            END DO
196!!gm BUG end
197         ENDIF
198                                                            !
199         
200      ELSE                                             !* Leap-Frog : Asselin filter and swap
201         !                                ! =============!
202         IF( ln_linssh ) THEN             ! Fixed volume !
203            !                             ! =============!
204            DO jk = 1, jpkm1                             
205               DO jj = 1, jpj
206                  DO ji = 1, jpi   
207                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
208                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
209                     !
210                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
211                     vb(ji,jj,jk) = zvf
212                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
213                     vn(ji,jj,jk) = va(ji,jj,jk)
214                  END DO
215               END DO
216            END DO
217            !                             ! ================!
218         ELSE                             ! Variable volume !
219            !                             ! ================!
220            ! Before scale factor at t-points
221            ! (used as a now filtered scale factor until the swap)
222            ! ----------------------------------------------------
223            DO jk = 1, jpkm1
224               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
225            END DO
226            ! Add volume filter correction: compatibility with tracer advection scheme
227            ! => time filter + conservation correction (only at the first level)
228            zcoef = atfp * rdt * r1_rau0
229
230            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
231
232            IF ( ln_rnf ) THEN
233               IF( ln_rnf_depth ) THEN
234                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
235                     DO jj = 1, jpj
236                        DO ji = 1, jpi
237                           IF( 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            END IF
248
249            IF ( ln_isf ) THEN   ! if ice shelf melting
250               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
251                  DO jj = 1, jpj
252                     DO ji = 1, jpi
253                        IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN
254                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
255                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
256                        ELSEIF ( jk==misfkb(ji,jj) ) THEN
257                           e3t_b(ji,jj,jk) = e3t_b(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
258                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
259                        ENDIF
260                     END DO
261                  END DO
262               END DO
263            END IF
264            !
265            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
266               ! Before filtered scale factor at (u/v)-points
267               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
268               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
269               DO jk = 1, jpkm1
270                  DO jj = 1, jpj
271                     DO ji = 1, jpi
272                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
273                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
274                        !
275                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
276                        vb(ji,jj,jk) = zvf
277                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
278                        vn(ji,jj,jk) = va(ji,jj,jk)
279                     END DO
280                  END DO
281               END DO
282               !
283            ELSE                          ! Asselin filter applied on thickness weighted velocity
284               !
285               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
286               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
287               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
288               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
289               DO jk = 1, jpkm1
290                  DO jj = 1, jpj
291                     DO ji = 1, jpi                 
292                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
293                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
294                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
295                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
296                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
297                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
298                        !
299                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
300                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
301                        !
302                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
303                        vb(ji,jj,jk) = zvf
304                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
305                        vn(ji,jj,jk) = va(ji,jj,jk)
306                     END DO
307                  END DO
308               END DO
309               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
310               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
311               !
312               DEALLOCATE( ze3u_f , ze3v_f )
313            ENDIF
314            !
315         ENDIF
316         !
317         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
318            ! Revert "before" velocities to time split estimate
319            ! Doing it here also means that asselin filter contribution is removed 
320            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
321            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
322            DO jk = 2, jpkm1
323               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
324               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
325            END DO
326            DO jk = 1, jpkm1
327               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
328               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
329            END DO
330         ENDIF
331         !
332      ENDIF ! neuler =/0
333      !
334      ! Set "now" and "before" barotropic velocities for next time step:
335      ! JC: Would be more clever to swap variables than to make a full vertical
336      ! integration
337      !
338      !
339      IF(.NOT.ln_linssh ) THEN
340         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
341         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
342         DO jk = 2, jpkm1
343            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
344            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
345         END DO
346         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
347         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
348      ENDIF
349      !
350      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
351      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
352      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
353      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
354      DO jk = 2, jpkm1
355         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
356         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
357         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
358         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
359      END DO
360      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
361      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
362      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
363      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
364      !
365      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
366         CALL iom_put(  "ubar", un_b(:,:) )
367         CALL iom_put(  "vbar", vn_b(:,:) )
368      ENDIF
369      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
370         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
371         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
372         CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm )
373      ENDIF
374      !
375      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
376         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
377      !
378      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
379      IF( l_trddyn     )   DEALLOCATE( zua, zva )
380      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
381      !
382   END SUBROUTINE dyn_nxt
383
384   !!=========================================================================
385END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.