source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 4 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 17.3 KB
RevLine 
[3]1MODULE dynnxt
[1502]2   !!=========================================================================
[3]3   !!                       ***  MODULE  dynnxt  ***
4   !! Ocean dynamics: time stepping
[1502]5   !!=========================================================================
[1438]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.
[1502]16   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
[2528]17   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
[2723]18   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
[4292]19   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
[6140]20   !!            3.6  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
[5930]21   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification
[1502]22   !!-------------------------------------------------------------------------
[1438]23 
[1502]24   !!-------------------------------------------------------------------------
[6140]25   !!   dyn_nxt       : obtain the next (after) horizontal velocity
[1502]26   !!-------------------------------------------------------------------------
[6140]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 dynadv         ! dynamics: vector invariant versus flux form
32   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
33   USE domvvl         ! variable volume
[7646]34   USE bdy_oce   , ONLY: ln_bdy
[6140]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
[4990]41   !
[6140]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 prtctl         ! Print control
47   USE timing         ! Timing
[2528]48#if defined key_agrif
49   USE agrif_opa_interp
50#endif
[3]51
52   IMPLICIT NONE
53   PRIVATE
54
[1438]55   PUBLIC    dyn_nxt   ! routine called by step.F90
56
[2715]57   !!----------------------------------------------------------------------
[2528]58   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1438]59   !! $Id$
[2715]60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
[3]62CONTAINS
63
64   SUBROUTINE dyn_nxt ( kt )
65      !!----------------------------------------------------------------------
66      !!                  ***  ROUTINE dyn_nxt  ***
67      !!                   
[5930]68      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
69      !!             condition on the after velocity, achieve the time stepping
[1502]70      !!             by applying the Asselin filter on now fields and swapping
71      !!             the fields.
[3]72      !!
[5930]73      !! ** Method  : * Ensure after velocities transport matches time splitting
74      !!             estimate (ln_dynspg_ts=T)
[3]75      !!
[1502]76      !!              * Apply lateral boundary conditions on after velocity
77      !!             at the local domain boundaries through lbc_lnk call,
[7646]78      !!             at the one-way open boundaries (ln_bdy=T),
[4990]79      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
[3]80      !!
[1502]81      !!              * Apply the time filter applied and swap of the dynamics
82      !!             arrays to start the next time step:
83      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
84      !!                (un,vn) = (ua,va).
[6140]85      !!             Note that with flux form advection and non linear free surface,
86      !!             the time filter is applied on thickness weighted velocity.
87      !!             As a result, dyn_nxt MUST be called after tra_nxt.
[1502]88      !!
89      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
90      !!               un,vn   now horizontal velocity of next time-step
[3]91      !!----------------------------------------------------------------------
92      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
[2715]93      !
[3]94      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[6140]95      INTEGER  ::   ikt          ! local integers
96      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
[4990]97      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
[7910]98      REAL(wp), DIMENSION(jpi,jpj)   ::  zue, zve
99      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ze3u_f, ze3v_f, zua, zva 
[1502]100      !!----------------------------------------------------------------------
[3294]101      !
[4990]102      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
[3294]103      !
104      !
[3]105      IF( kt == nit000 ) THEN
106         IF(lwp) WRITE(numout,*)
107         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
108         IF(lwp) WRITE(numout,*) '~~~~~~~'
109      ENDIF
110
[5930]111      IF ( ln_dynspg_ts ) THEN
112         ! Ensure below that barotropic velocities match time splitting estimate
113         ! Compute actual transport and replace it with ts estimate at "after" time step
[7753]114         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
115         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
[5930]116         DO jk = 2, jpkm1
[7753]117            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
118            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
[1502]119         END DO
120         DO jk = 1, jpkm1
[7753]121            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
122            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
[592]123         END DO
[6140]124         !
125         IF( .NOT.ln_bt_fw ) THEN
[5930]126            ! Remove advective velocity from "now velocities"
127            ! prior to asselin filtering     
128            ! In the forward case, this is done below after asselin filtering   
129            ! so that asselin contribution is removed at the same time
130            DO jk = 1, jpkm1
[7753]131               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk)
132               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk)
133            END DO 
[5930]134         ENDIF
[4292]135      ENDIF
136
[1502]137      ! Update after velocity on domain lateral boundaries
138      ! --------------------------------------------------     
[5930]139# if defined key_agrif
140      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
141# endif
142      !
[1502]143      CALL lbc_lnk( ua, 'U', -1. )     !* local domain boundaries
144      CALL lbc_lnk( va, 'V', -1. ) 
145      !
146      !                                !* BDY open boundaries
[7646]147      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
148      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
[3294]149
150!!$   Do we need a call to bdy_vol here??
151      !
[4990]152      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
153         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
154         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
155         !
156         !                                  ! Kinetic energy and Conversion
157         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
158         !
159         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
[7753]160            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
161            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
[4990]162            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
163            CALL iom_put( "vtrd_tot", zva )
164         ENDIF
165         !
[7753]166         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
167         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
168         !                                  !  computation of the asselin filter trends)
[4990]169      ENDIF
170
[1438]171      ! Time filter and swap of dynamics arrays
172      ! ------------------------------------------
[1502]173      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
174         DO jk = 1, jpkm1
[7753]175            un(:,:,jk) = ua(:,:,jk)                          ! un <-- ua
176            vn(:,:,jk) = va(:,:,jk)
[1438]177         END DO
[6140]178         IF(.NOT.ln_linssh ) THEN
[4292]179            DO jk = 1, jpkm1
[7753]180               e3t_b(:,:,jk) = e3t_n(:,:,jk)
181               e3u_b(:,:,jk) = e3u_n(:,:,jk)
182               e3v_b(:,:,jk) = e3v_n(:,:,jk)
[6140]183            END DO
[4292]184         ENDIF
[1502]185      ELSE                                             !* Leap-Frog : Asselin filter and swap
[2528]186         !                                ! =============!
[6140]187         IF( ln_linssh ) THEN             ! Fixed volume !
[2528]188            !                             ! =============!
[1502]189            DO jk = 1, jpkm1                             
[592]190               DO jj = 1, jpj
[1502]191                  DO ji = 1, jpi   
[4990]192                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
193                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
[1502]194                     !
195                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
196                     vb(ji,jj,jk) = zvf
197                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
198                     vn(ji,jj,jk) = va(ji,jj,jk)
199                  END DO
200               END DO
201            END DO
[2528]202            !                             ! ================!
203         ELSE                             ! Variable volume !
204            !                             ! ================!
[4292]205            ! Before scale factor at t-points
206            ! (used as a now filtered scale factor until the swap)
207            ! ----------------------------------------------------
[6140]208            IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN    ! No asselin filtering on thicknesses if forward time splitting
[7753]209               e3t_b(:,:,1:jpkm1) = e3t_n(:,:,1:jpkm1)
[4292]210            ELSE
[6140]211               DO jk = 1, jpkm1
[7753]212                  e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
[6140]213               END DO
[4292]214               ! Add volume filter correction: compatibility with tracer advection scheme
215               ! => time filter + conservation correction (only at the first level)
[6140]216               zcoef = atfp * rdt * r1_rau0
217               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting
[7753]218                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) &
219                                 &                      - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)
[5643]220               ELSE                     ! if ice shelf melting
[6140]221                  DO jj = 1, jpj
222                     DO ji = 1, jpi
223                        ikt = mikt(ji,jj)
224                        e3t_b(ji,jj,ikt) = e3t_b(ji,jj,ikt) - zcoef * (  emp_b   (ji,jj) - emp   (ji,jj)  &
225                           &                                           - rnf_b   (ji,jj) + rnf   (ji,jj)  &
226                           &                                           + fwfisf_b(ji,jj) - fwfisf(ji,jj)  ) * tmask(ji,jj,ikt)
[5643]227                     END DO
228                  END DO
229               END IF
[4292]230            ENDIF
[2528]231            !
[6140]232            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
233               ! Before filtered scale factor at (u/v)-points
234               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
235               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
[4292]236               DO jk = 1, jpkm1
237                  DO jj = 1, jpj
[2528]238                     DO ji = 1, jpi
[4292]239                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
240                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
[2528]241                        !
242                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
243                        vb(ji,jj,jk) = zvf
244                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
245                        vn(ji,jj,jk) = va(ji,jj,jk)
246                     END DO
247                  END DO
248               END DO
249               !
[6140]250            ELSE                          ! Asselin filter applied on thickness weighted velocity
251               !
252               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
253               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
254               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
[4292]255               DO jk = 1, jpkm1
256                  DO jj = 1, jpj
[4312]257                     DO ji = 1, jpi                 
[6140]258                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
259                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
260                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
261                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
262                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
263                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
[2528]264                        !
[3294]265                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
266                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
[2528]267                        !
[3294]268                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
[2528]269                        vb(ji,jj,jk) = zvf
[3294]270                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
[2528]271                        vn(ji,jj,jk) = va(ji,jj,jk)
272                     END DO
273                  END DO
274               END DO
[7753]275               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
276               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
[6140]277               !
[2528]278            ENDIF
279            !
[3]280         ENDIF
[2528]281         !
[6140]282         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
[4312]283            ! Revert "before" velocities to time split estimate
284            ! Doing it here also means that asselin filter contribution is removed 
[7753]285            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
286            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)   
[4990]287            DO jk = 2, jpkm1
[7753]288               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
289               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)   
[4370]290            END DO
291            DO jk = 1, jpkm1
[7753]292               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
293               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
[4292]294            END DO
295         ENDIF
296         !
297      ENDIF ! neuler =/0
[4354]298      !
299      ! Set "now" and "before" barotropic velocities for next time step:
300      ! JC: Would be more clever to swap variables than to make a full vertical
301      ! integration
302      !
[4370]303      !
[6140]304      IF(.NOT.ln_linssh ) THEN
[7753]305         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
306         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
[6140]307         DO jk = 2, jpkm1
[7753]308            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
309            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
[4354]310         END DO
[7753]311         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
312         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
[4354]313      ENDIF
314      !
[7753]315      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
316      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
317      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
318      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
[6140]319      DO jk = 2, jpkm1
[7753]320         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
321         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
322         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
323         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
[4354]324      END DO
[7753]325      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
326      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
327      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
328      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
[4354]329      !
[6140]330      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
331         CALL iom_put(  "ubar", un_b(:,:) )
332         CALL iom_put(  "vbar", vn_b(:,:) )
333      ENDIF
[4990]334      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
[7753]335         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
336         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
[4990]337         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
338      ENDIF
339      !
[1438]340      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
341         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
[6140]342      !
[2715]343      !
[3294]344      IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
345      !
[3]346   END SUBROUTINE dyn_nxt
347
[1502]348   !!=========================================================================
[3]349END MODULE dynnxt
Note: See TracBrowser for help on using the repository browser.