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

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

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps: tidy up.
  • Property svn:keywords set to Id
File size: 18.4 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, kNm1, kNp1, puu, pvv, pe3t, pe3u, pe3v )
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 ) ::   kNm1, kNp1 ! before and after time level indices
96      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: puu, pvv ! now velocities to be time filtered
97      REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pe3t, pe3u, pe3v ! now scale factors to be time filtered
98      !
99      INTEGER  ::   ji, jj, jk   ! dummy loop indices
100      INTEGER  ::   ikt          ! local integers
101      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars
102      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      -
103      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
104      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva 
105      !!----------------------------------------------------------------------
106      !
107      IF( ln_timing    )   CALL timing_start('dyn_nxt')
108      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
109      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
110      !
111      IF( kt == nit000 ) THEN
112         IF(lwp) WRITE(numout,*)
113         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
114         IF(lwp) WRITE(numout,*) '~~~~~~~'
115      ENDIF
116
117      IF ( ln_dynspg_ts ) THEN
118         ! Ensure below that barotropic velocities match time splitting estimate
119         ! Compute actual transport and replace it with ts estimate at "after" time step
120         zue(:,:) = e3u(:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1)
121         zve(:,:) = e3v(:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1)
122         DO jk = 2, jpkm1
123            zue(:,:) = zue(:,:) + e3u(:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk)
124            zve(:,:) = zve(:,:) + e3v(:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk)
125         END DO
126         DO jk = 1, jpkm1
127            uu(:,:,jk,kNp1) = ( uu(:,:,jk,kNp1) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
128            vv(:,:,jk,kNp1) = ( vv(:,:,jk,kNp1) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
129         END DO
130         !
131         IF( .NOT.ln_bt_fw ) THEN
132            ! Remove advective velocity from "now velocities"
133            ! prior to asselin filtering     
134            ! In the forward case, this is done below after asselin filtering   
135            ! so that asselin contribution is removed at the same time
136            DO jk = 1, jpkm1
137               puu(:,:,jk) = ( puu(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
138               pvv(:,:,jk) = ( pvv(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
139            END DO 
140         ENDIF
141      ENDIF
142
143      ! Update after velocity on domain lateral boundaries
144      ! --------------------------------------------------     
145# if defined key_agrif
146      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
147# endif
148      !
149      CALL lbc_lnk_multi( 'dynnxt', uu(:,:,:,kNp1), 'U', -1., vv(:,:,:,kNp1), 'V', -1. )     !* local domain boundaries
150      !
151      !                                !* BDY open boundaries
152      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
153      IF( ln_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      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
158         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
159         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
160         !
161         !                                  ! Kinetic energy and Conversion
162         IF( ln_KE_trd  )   CALL trd_dyn( uu(:,:,:,kNp1), vv(:,:,:,kNp1), jpdyn_ken, kt )
163         !
164         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
165            zua(:,:,:) = ( uu(:,:,:,kNp1) - uu(:,:,:,kNm1) ) * z1_2dt
166            zva(:,:,:) = ( vv(:,:,:,kNp1) - vv(:,:,:,kNm1) ) * z1_2dt
167            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
168            CALL iom_put( "vtrd_tot", zva )
169         ENDIF
170         !
171         zua(:,:,:) = puu(:,:,:)             ! save the now velocity before the asselin filter
172         zva(:,:,:) = pvv(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
173         !                                  !  computation of the asselin filter trends)
174      ENDIF
175
176      ! Time filter and swap of dynamics arrays
177      ! ------------------------------------------
178         
179      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin filter and swap
180         !                                ! =============!
181         IF( ln_linssh ) THEN             ! Fixed volume !
182            !                             ! =============!
183            DO jk = 1, jpkm1                             
184               DO jj = 1, jpj
185                  DO ji = 1, jpi   
186                     puu(ji,jj,jk) = puu(ji,jj,jk) + atfp * ( uu(ji,jj,jk,kNm1) - 2._wp * puu(ji,jj,jk) + uu(ji,jj,jk,kNp1) )
187                     pvv(ji,jj,jk) = pvv(ji,jj,jk) + atfp * ( vv(ji,jj,jk,kNm1) - 2._wp * pvv(ji,jj,jk) + vv(ji,jj,jk,kNp1) )
188                  END DO
189               END DO
190            END DO
191            !                             ! ================!
192         ELSE                             ! Variable volume !
193            !                             ! ================!
194            ! Time-filtered scale factor at t-points
195            ! ----------------------------------------------------
196            ALLOCATE( ze3t_f(jpi,jpj,jpk) )
197            DO jk = 1, jpkm1
198               ze3t_f(:,:,jk) = pe3t(:,:,jk) + atfp * ( e3t(:,:,jk,kNm1) - 2._wp * pe3t(:,:,jk) + e3t(:,:,jk,kNp1) )
199            END DO
200            ! Add volume filter correction: compatibility with tracer advection scheme
201            ! => time filter + conservation correction (only at the first level)
202            zcoef = atfp * rdt * r1_rau0
203
204            ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
205
206            IF ( ln_rnf ) THEN
207               IF( ln_rnf_depth ) THEN
208                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
209                     DO jj = 1, jpj
210                        DO ji = 1, jpi
211                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
212                               ze3t_f(ji,jj,jk) =   ze3t_f(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
213                                      &          * ( pe3t(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
214                           ENDIF
215                        ENDDO
216                     ENDDO
217                  ENDDO
218               ELSE
219                  ze3t_f(:,:,1) = ze3t_f(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
220               ENDIF
221            END IF
222
223            IF ( ln_isf ) THEN   ! if ice shelf melting
224               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
225                  DO jj = 1, jpj
226                     DO ji = 1, jpi
227                        IF( misfkt(ji,jj) <=jk .and. jk < misfkb(ji,jj)  ) THEN
228                           ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
229                                &          * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
230                        ELSEIF ( jk==misfkb(ji,jj) ) THEN
231                           ze3t_f(ji,jj,jk) = ze3t_f(ji,jj,jk) - zcoef * ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
232                                &          * ( pe3t(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * ralpha(ji,jj) * tmask(ji,jj,jk)
233                        ENDIF
234                     END DO
235                  END DO
236               END DO
237            END IF
238            !
239            pe3t(:,:,1:jpkm1) = ze3t_f(:,:,1:jpkm1)        ! filtered scale factor at T-points
240            !
241            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
242               ! Before filtered scale factor at (u/v)-points
243               CALL dom_vvl_interpol( pe3t(:,:,:), pe3u(:,:,:), 'U' )
244               CALL dom_vvl_interpol( pe3t(:,:,:), pe3v(:,:,:), 'V' )
245               DO jk = 1, jpkm1
246                  DO jj = 1, jpj
247                     DO ji = 1, jpi
248                        puu(ji,jj,jk) = puu(ji,jj,jk) + atfp * ( uu(ji,jj,jk,kNm1) - 2._wp * puu(ji,jj,jk) + uu(ji,jj,jk,kNp1) )
249                        pvv(ji,jj,jk) = pvv(ji,jj,jk) + atfp * ( vv(ji,jj,jk,kNm1) - 2._wp * pvv(ji,jj,jk) + vv(ji,jj,jk,kNp1) )
250                     END DO
251                  END DO
252               END DO
253               !
254            ELSE                          ! Asselin filter applied on thickness weighted velocity
255               !
256               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
257               ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
258               CALL dom_vvl_interpol( pe3t(:,:,:), ze3u_f, 'U' )
259               CALL dom_vvl_interpol( pe3t(:,:,:), ze3v_f, 'V' )
260               DO jk = 1, jpkm1
261                  DO jj = 1, jpj
262                     DO ji = 1, jpi                 
263                        zue3a = e3u(ji,jj,jk,kNp1) * uu(ji,jj,jk,kNp1)
264                        zve3a = e3v(ji,jj,jk,kNp1) * vv(ji,jj,jk,kNp1)
265                        zue3n = pe3u(ji,jj,jk)     * puu(ji,jj,jk)
266                        zve3n = pe3v(ji,jj,jk)     * pvv(ji,jj,jk)
267                        zue3b = e3u(ji,jj,jk,kNm1) * uu(ji,jj,jk,kNm1)
268                        zve3b = e3v(ji,jj,jk,kNm1) * vv(ji,jj,jk,kNm1)
269                        !
270                        puu(ji,jj,jk) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
271                        pvv(ji,jj,jk) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
272                     END DO
273                  END DO
274               END DO
275               pe3u(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) 
276               pe3v(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
277               !
278               DEALLOCATE( ze3u_f , ze3v_f )
279            ENDIF
280            !
281            DEALLOCATE( ze3t_f )
282         ENDIF
283         !
284         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
285            ! Revert filtered "now" velocities to time split estimate
286            ! Doing it here also means that asselin filter contribution is removed 
287            zue(:,:) = pe3u(:,:,1) * puu(:,:,1) * umask(:,:,1)
288            zve(:,:) = pe3v(:,:,1) * pvv(:,:,1) * vmask(:,:,1)   
289            DO jk = 2, jpkm1
290               zue(:,:) = zue(:,:) + pe3u(:,:,jk) * puu(:,:,jk) * umask(:,:,jk)
291               zve(:,:) = zve(:,:) + pe3v(:,:,jk) * pvv(:,:,jk) * vmask(:,:,jk)   
292            END DO
293            DO jk = 1, jpkm1
294               puu(:,:,jk) = puu(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
295               pvv(:,:,jk) = pvv(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
296            END DO
297         ENDIF
298         !
299         ikt = Nnn
300      ELSE
301         ikt = kNm1
302         puu(:,:,:)  =  uu(:,:,:,kNp1)
303         pvv(:,:,:)  =  vv(:,:,:,kNp1)
304         pe3t(:,:,:) = e3t(:,:,:,kNp1)
305         pe3u(:,:,:) = e3u(:,:,:,kNp1)
306         pe3v(:,:,:) = e3v(:,:,:,kNp1)
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      ! DS IMMERSE :: This is very confusing as it stands. But should
314      ! hopefully be simpler when we do the time-level swapping for the
315      ! external mode solution.
316      !
317      IF(.NOT.ln_linssh ) THEN
318         hu(:,:,ikt) = e3u(:,:,1,ikt ) * umask(:,:,1)
319         hv(:,:,ikt) = e3v(:,:,1,ikt ) * vmask(:,:,1)
320         DO jk = 2, jpkm1
321            hu(:,:,ikt) = hu(:,:,ikt) + e3u(:,:,jk,ikt ) * umask(:,:,jk)
322            hv(:,:,ikt) = hv(:,:,ikt) + e3v(:,:,jk,ikt ) * vmask(:,:,jk)
323         END DO
324         r1_hu(:,:,ikt) = ssumask(:,:) / ( hu(:,:,ikt) + 1._wp - ssumask(:,:) )
325         r1_hv(:,:,ikt) = ssvmask(:,:) / ( hv(:,:,ikt) + 1._wp - ssvmask(:,:) )
326      ENDIF
327      !
328      un_b(:,:) = e3u(:,:,1,kNp1) * uu(:,:,1,kNp1) * umask(:,:,1)
329      ub_b(:,:) = e3u(:,:,1,ikt ) * uu(:,:,1,ikt ) * umask(:,:,1)
330      vn_b(:,:) = e3v(:,:,1,kNp1) * vv(:,:,1,kNp1) * vmask(:,:,1)
331      vb_b(:,:) = e3v(:,:,1,ikt ) * vv(:,:,1,ikt ) * vmask(:,:,1)
332      DO jk = 2, jpkm1
333         un_b(:,:) = un_b(:,:) + e3u(:,:,jk,kNp1) * uu(:,:,jk,kNp1) * umask(:,:,jk)
334         ub_b(:,:) = ub_b(:,:) + e3u(:,:,jk,ikt ) * uu(:,:,jk,ikt ) * umask(:,:,jk)
335         vn_b(:,:) = vn_b(:,:) + e3v(:,:,jk,kNp1) * vv(:,:,jk,kNp1) * vmask(:,:,jk)
336         vb_b(:,:) = vb_b(:,:) + e3v(:,:,jk,ikt ) * vv(:,:,jk,ikt ) * vmask(:,:,jk)
337      END DO
338      un_b(:,:) = un_b(:,:) * r1_hu(:,:,kNp1)
339      vn_b(:,:) = vn_b(:,:) * r1_hv(:,:,kNp1)
340      ub_b(:,:) = ub_b(:,:) * r1_hu(:,:,ikt)
341      vb_b(:,:) = vb_b(:,:) * r1_hv(:,:,ikt)
342      !
343      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
344         CALL iom_put(  "ubar", un_b(:,:) )
345         CALL iom_put(  "vbar", vn_b(:,:) )
346      ENDIF
347      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
348         zua(:,:,:) = ( puu(:,:,:) - zua(:,:,:) ) * z1_2dt
349         zva(:,:,:) = ( pvv(:,:,:) - zva(:,:,:) ) * z1_2dt
350         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
351      ENDIF
352      !
353      IF(ln_ctl)   CALL prt_ctl( tab3d_1=uu(:,:,:,kNp1), clinfo1=' nxt  - uu(:,:,:,kNp1): ', mask1=umask,   &
354         &                       tab3d_2=vv(:,:,:,kNp1), clinfo2=' vv(:,:,:,kNp1): '       , mask2=vmask )
355      !
356      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
357      IF( l_trddyn     )   DEALLOCATE( zua, zva )
358      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
359      !
360   END SUBROUTINE dyn_nxt
361
362   !!=========================================================================
363END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.