source: NEMO/trunk/src/OCE/DYN/dynatf.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 9 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

  • Property svn:keywords set to Id
File size: 17.2 KB
Line 
1MODULE dynatf
2   !!=========================================================================
3   !!                       ***  MODULE  dynatf  ***
4   !! Ocean dynamics: time filtering
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   !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf.F90. Now just does time filtering.
23   !!-------------------------------------------------------------------------
24 
25   !!----------------------------------------------------------------------------------------------
26   !!   dyn_atf       : apply Asselin time filtering to "now" velocities and vertical scale factors
27   !!----------------------------------------------------------------------------------------------
28   USE oce            ! ocean dynamics and tracers
29   USE dom_oce        ! ocean space and time domain
30   USE sbc_oce        ! Surface boundary condition: ocean fields
31   USE sbcrnf         ! river runoffs
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   USE isf_oce   , ONLY: ln_isf     ! ice shelf
44   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine
45   !
46   USE in_out_manager ! I/O manager
47   USE iom            ! I/O manager library
48   USE lbclnk         ! lateral boundary condition (or mpp link)
49   USE lib_mpp        ! MPP library
50   USE prtctl         ! Print control
51   USE timing         ! Timing
52#if defined key_agrif
53   USE agrif_oce_interp
54#endif
55
56   IMPLICIT NONE
57   PRIVATE
58
59   PUBLIC    dyn_atf   ! routine called by step.F90
60
61   !! * Substitutions
62#  include "do_loop_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
65   !! $Id$
66   !! Software governed by the CeCILL license (see ./LICENSE)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v )
71      !!----------------------------------------------------------------------
72      !!                  ***  ROUTINE dyn_atf  ***
73      !!                   
74      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
75      !!             condition on the after velocity and apply the Asselin time
76      !!             filter to the now fields.
77      !!
78      !! ** Method  : * Ensure after velocities transport matches time splitting
79      !!             estimate (ln_dynspg_ts=T)
80      !!
81      !!              * Apply lateral boundary conditions on after velocity
82      !!             at the local domain boundaries through lbc_lnk call,
83      !!             at the one-way open boundaries (ln_bdy=T),
84      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
85      !!
86      !!              * Apply the Asselin time filter to the now fields
87      !!             arrays to start the next time step:
88      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm))
89      !!                                    + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ]
90      !!             Note that with flux form advection and non linear free surface,
91      !!             the time filter is applied on thickness weighted velocity.
92      !!             As a result, dyn_atf MUST be called after tra_atf.
93      !!
94      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity
95      !!----------------------------------------------------------------------
96      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
97      INTEGER                             , INTENT(in   ) :: Kbb, Kmm, Kaa    ! before and after time level indices
98      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered
99      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered
100      !
101      INTEGER  ::   ji, jj, jk   ! dummy loop indices
102      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars
103      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      -
104      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve, zwfld
105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva 
106      !!----------------------------------------------------------------------
107      !
108      IF( ln_timing    )   CALL timing_start('dyn_atf')
109      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
110      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
111      !
112      IF( kt == nit000 ) THEN
113         IF(lwp) WRITE(numout,*)
114         IF(lwp) WRITE(numout,*) 'dyn_atf : Asselin time filtering'
115         IF(lwp) WRITE(numout,*) '~~~~~~~'
116      ENDIF
117
118      IF ( ln_dynspg_ts ) THEN
119         ! Ensure below that barotropic velocities match time splitting estimate
120         ! Compute actual transport and replace it with ts estimate at "after" time step
121         zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
122         zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
123         DO jk = 2, jpkm1
124            zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
125            zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
126         END DO
127         DO jk = 1, jpkm1
128            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
129            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
130         END DO
131         !
132         IF( .NOT.ln_bt_fw ) THEN
133            ! Remove advective velocity from "now velocities"
134            ! prior to asselin filtering     
135            ! In the forward case, this is done below after asselin filtering   
136            ! so that asselin contribution is removed at the same time
137            DO jk = 1, jpkm1
138               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
139               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
140            END DO 
141         ENDIF
142      ENDIF
143
144      ! Update after velocity on domain lateral boundaries
145      ! --------------------------------------------------     
146# if defined key_agrif
147      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
148# endif
149      !
150      CALL lbc_lnk_multi( 'dynatf', puu(:,:,:,Kaa), 'U', -1., pvv(:,:,:,Kaa), 'V', -1. )     !* local domain boundaries
151      !
152      !                                !* BDY open boundaries
153      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
154      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
155
156!!$   Do we need a call to bdy_vol here??
157      !
158      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
159         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
160         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
161         !
162         !                                  ! Kinetic energy and Conversion
163         IF( ln_KE_trd  )   CALL trd_dyn( puu(:,:,:,Kaa), pvv(:,:,:,Kaa), jpdyn_ken, kt, Kmm )
164         !
165         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
166            zua(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) * z1_2dt
167            zva(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) * z1_2dt
168            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
169            CALL iom_put( "vtrd_tot", zva )
170         ENDIF
171         !
172         zua(:,:,:) = puu(:,:,:,Kmm)             ! save the now velocity before the asselin filter
173         zva(:,:,:) = pvv(:,:,:,Kmm)             ! (caution: there will be a shift by 1 timestep in the
174         !                                  !  computation of the asselin filter trends)
175      ENDIF
176
177      ! Time filter and swap of dynamics arrays
178      ! ------------------------------------------
179         
180      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter
181         !                                ! =============!
182         IF( ln_linssh ) THEN             ! Fixed volume !
183            !                             ! =============!
184            DO_3D_11_11( 1, jpkm1 )
185               puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
186               pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
187            END_3D
188            !                             ! ================!
189         ELSE                             ! Variable volume !
190            !                             ! ================!
191            ! Time-filtered scale factor at t-points
192            ! ----------------------------------------------------
193            ALLOCATE( ze3t_f(jpi,jpj,jpk), zwfld(jpi,jpj) )
194            DO jk = 1, jpkm1
195               ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) )
196            END DO
197            ! Add volume filter correction: compatibility with tracer advection scheme
198            ! => time filter + conservation correction
199            zcoef = atfp * rdt * r1_rau0
200            zwfld(:,:) = emp_b(:,:) - emp(:,:)
201            IF ( ln_rnf ) zwfld(:,:) =  zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) )
202            DO jk = 1, jpkm1
203               ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) &
204                              &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 
205            END DO
206            !
207            ! ice shelf melting (deal separately as it can be in depth)
208            ! PM: we could probably define a generic subroutine to do the in depth correction
209            !     to manage rnf, isf and possibly in the futur icb, tide water glacier (...)
210            !     ...(kt, coef, ktop, kbot, hz, fwf_b, fwf)
211            IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, atfp * rdt )
212            !
213            pe3t(:,:,1:jpkm1,Kmm) = ze3t_f(:,:,1:jpkm1)        ! filtered scale factor at T-points
214            !
215            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
216               ! Before filtered scale factor at (u/v)-points
217               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' )
218               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' )
219               DO_3D_11_11( 1, jpkm1 )
220                  puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
221                  pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
222               END_3D
223               !
224            ELSE                          ! Asselin filter applied on thickness weighted velocity
225               !
226               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
227               ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
228               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' )
229               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' )
230               DO_3D_11_11( 1, jpkm1 )
231                  zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa)
232                  zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa)
233                  zue3n = pe3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
234                  zve3n = pe3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
235                  zue3b = pe3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb)
236                  zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb)
237                  !
238                  puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
239                  pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
240               END_3D
241               pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 
242               pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1)
243               !
244               DEALLOCATE( ze3u_f , ze3v_f )
245            ENDIF
246            !
247            DEALLOCATE( ze3t_f, zwfld )
248         ENDIF
249         !
250         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
251            ! Revert filtered "now" velocities to time split estimate
252            ! Doing it here also means that asselin filter contribution is removed 
253            zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
254            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)   
255            DO jk = 2, jpkm1
256               zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
257               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)   
258            END DO
259            DO jk = 1, jpkm1
260               puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) - (zue(:,:) * r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)) * umask(:,:,jk)
261               pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) - (zve(:,:) * r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)) * vmask(:,:,jk)
262            END DO
263         ENDIF
264         !
265      ENDIF ! neuler /= 0
266      !
267      ! Set "now" and "before" barotropic velocities for next time step:
268      ! JC: Would be more clever to swap variables than to make a full vertical
269      ! integration
270      !
271      IF(.NOT.ln_linssh ) THEN
272         hu(:,:,Kmm) = pe3u(:,:,1,Kmm ) * umask(:,:,1)
273         hv(:,:,Kmm) = pe3v(:,:,1,Kmm ) * vmask(:,:,1)
274         DO jk = 2, jpkm1
275            hu(:,:,Kmm) = hu(:,:,Kmm) + pe3u(:,:,jk,Kmm ) * umask(:,:,jk)
276            hv(:,:,Kmm) = hv(:,:,Kmm) + pe3v(:,:,jk,Kmm ) * vmask(:,:,jk)
277         END DO
278         r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
279         r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
280      ENDIF
281      !
282      uu_b(:,:,Kaa) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
283      uu_b(:,:,Kmm) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
284      vv_b(:,:,Kaa) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
285      vv_b(:,:,Kmm) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
286      DO jk = 2, jpkm1
287         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
288         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
289         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
290         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
291      END DO
292      uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa)
293      vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa)
294      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
295      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
296      !
297      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
298         CALL iom_put(  "ubar", uu_b(:,:,Kmm) )
299         CALL iom_put(  "vbar", vv_b(:,:,Kmm) )
300      ENDIF
301      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
302         zua(:,:,:) = ( puu(:,:,:,Kmm) - zua(:,:,:) ) * z1_2dt
303         zva(:,:,:) = ( pvv(:,:,:,Kmm) - zva(:,:,:) ) * z1_2dt
304         CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm )
305      ENDIF
306      !
307      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   &
308         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask )
309      !
310      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
311      IF( l_trddyn     )   DEALLOCATE( zua, zva )
312      IF( ln_timing    )   CALL timing_stop('dyn_atf')
313      !
314   END SUBROUTINE dyn_atf
315
316   !!=========================================================================
317END MODULE dynatf
Note: See TracBrowser for help on using the repository browser.