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.
dynatfLF.F90 in NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynatfLF.F90 @ 12590

Last change on this file since 12590 was 12583, checked in by techene, 4 years ago

OCE/DOM/domqe.F90: add gdep at time level Kbb in dom_qe_sf_update, OCE/DOM/domzgr_substitute.h90: create the substitute module, OCE/DYN/dynatfLF.F90, OCE/TRA/traatfLF.F90: move boundary condition management and agrif management from atf modules to OCE/steplf.F90, OCE/SBC/sbcice_cice.F90, ICE/iceistate.F90 : remove dom_vvl_interpol and replace by dom_vvl_zgr ?
  • Property svn:keywords set to Id
File size: 16.1 KB
RevLine 
[12581]1MODULE dynatfLF
[1502]2   !!=========================================================================
[12581]3   !!                       ***  MODULE  dynatfLF  ***
[11050]4   !! Ocean dynamics: time filtering
[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
[12581]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
[12581]22   !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatfLF.F90. Now just does time filtering.
[1502]23   !!-------------------------------------------------------------------------
[12581]24
[11050]25   !!----------------------------------------------------------------------------------------------
[12581]26   !!   dyn_atfLF       : apply Asselin time filtering to "now" velocities and vertical scale factors
[11050]27   !!----------------------------------------------------------------------------------------------
[6140]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
[9023]31   USE sbcrnf         ! river runoffs
[6140]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
[7646]36   USE bdy_oce   , ONLY: ln_bdy
[6140]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
[12150]43   USE isf_oce   , ONLY: ln_isf     ! ice shelf
[12581]44   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine
[4990]45   !
[6140]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
[2528]52#if defined key_agrif
[9570]53   USE agrif_oce_interp
[11050]54#endif
[3]55
56   IMPLICIT NONE
57   PRIVATE
58
[12581]59   PUBLIC    dyn_atf_lf   ! routine called by step.F90
[1438]60
[12340]61   !! * Substitutions
62#  include "do_loop_substitute.h90"
[2715]63   !!----------------------------------------------------------------------
[9598]64   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[12581]65   !! $Id$
[10068]66   !! Software governed by the CeCILL license (see ./LICENSE)
[2715]67   !!----------------------------------------------------------------------
[3]68CONTAINS
69
[12581]70   SUBROUTINE dyn_atf_lf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v )
[3]71      !!----------------------------------------------------------------------
[12581]72      !!                  ***  ROUTINE dyn_atf_lf  ***
73      !!
74      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
[11475]75      !!             condition on the after velocity and apply the Asselin time
76      !!             filter to the now fields.
[3]77      !!
[5930]78      !! ** Method  : * Ensure after velocities transport matches time splitting
79      !!             estimate (ln_dynspg_ts=T)
[3]80      !!
[12581]81      !!              * Apply lateral boundary conditions on after velocity
[1502]82      !!             at the local domain boundaries through lbc_lnk call,
[7646]83      !!             at the one-way open boundaries (ln_bdy=T),
[4990]84      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
[3]85      !!
[11475]86      !!              * Apply the Asselin time filter to the now fields
[1502]87      !!             arrays to start the next time step:
[12581]88      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm))
[11475]89      !!                                    + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ]
[6140]90      !!             Note that with flux form advection and non linear free surface,
91      !!             the time filter is applied on thickness weighted velocity.
[12581]92      !!             As a result, dyn_atf_lf MUST be called after tra_atf.
[1502]93      !!
[12581]94      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity
[3]95      !!----------------------------------------------------------------------
[11050]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
[2715]100      !
[3]101      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[11050]102      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars
103      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      -
[12581]104      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zua, zva
[1502]106      !!----------------------------------------------------------------------
[3294]107      !
[12581]108      IF( ln_timing    )   CALL timing_start('dyn_atf_lf')
[9019]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) )
[3294]111      !
[3]112      IF( kt == nit000 ) THEN
113         IF(lwp) WRITE(numout,*)
[12581]114         IF(lwp) WRITE(numout,*) 'dyn_atf_lf : Asselin time filtering'
[3]115         IF(lwp) WRITE(numout,*) '~~~~~~~'
116      ENDIF
117
[12583]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( 'dynatfLF', 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. )
[4292]155
[3294]156!!$   Do we need a call to bdy_vol here??
157      !
[4990]158      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
[12581]159         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
[4990]160         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
161         !
162         !                                  ! Kinetic energy and Conversion
[11050]163         IF( ln_KE_trd  )   CALL trd_dyn( puu(:,:,:,Kaa), pvv(:,:,:,Kaa), jpdyn_ken, kt, Kmm )
[4990]164         !
165         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
[11050]166            zua(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) * z1_2dt
167            zva(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) * z1_2dt
[4990]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         !
[11050]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
[7753]174         !                                  !  computation of the asselin filter trends)
[4990]175      ENDIF
176
[1438]177      ! Time filter and swap of dynamics arrays
178      ! ------------------------------------------
[12581]179
180      IF( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN    !* Leap-Frog : Asselin time filter
[2528]181         !                                ! =============!
[6140]182         IF( ln_linssh ) THEN             ! Fixed volume !
[2528]183            !                             ! =============!
[12340]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
[2528]188            !                             ! ================!
189         ELSE                             ! Variable volume !
190            !                             ! ================!
[11050]191            ! Time-filtered scale factor at t-points
[4292]192            ! ----------------------------------------------------
[12581]193            DO jk = 1, jpk                                          ! filtered scale factor at T-points
194               pe3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + r3t_f(:,:) * tmask(:,:,jk) )
[9023]195            END DO
[2528]196            !
[12150]197            !
[6140]198            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
199               ! Before filtered scale factor at (u/v)-points
[12581]200               DO jk = 1, jpk
201                  pe3u(:,:,jk,Kmm) = e3u_0(:,:,jk) * ( 1._wp + r3u_f(:,:) * umask(:,:,jk) )
202                  pe3v(:,:,jk,Kmm) = e3v_0(:,:,jk) * ( 1._wp + r3v_f(:,:) * vmask(:,:,jk) )
203               END DO
204               !
[12340]205               DO_3D_11_11( 1, jpkm1 )
206                  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) )
207                  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) )
208               END_3D
[2528]209               !
[6140]210            ELSE                          ! Asselin filter applied on thickness weighted velocity
211               !
[12340]212               DO_3D_11_11( 1, jpkm1 )
213                  zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa)
214                  zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa)
215                  zue3n = pe3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
216                  zve3n = pe3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
217                  zue3b = pe3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb)
218                  zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb)
[12581]219                  !                                                 ! filtered scale factor at U-,V-points
220                  pe3u(ji,jj,jk,Kmm) = e3u_0(ji,jj,jk) * ( 1._wp + r3u_f(ji,jj) * umask(ji,jj,jk) )
221                  pe3v(ji,jj,jk,Kmm) = e3v_0(ji,jj,jk) * ( 1._wp + r3v_f(ji,jj) * vmask(ji,jj,jk) )
[12340]222                  !
[12581]223                  puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / pe3u(ji,jj,jk,Kmm) !!st ze3u_f(ji,jj,jk)
224                  pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / pe3v(ji,jj,jk,Kmm) !!st ze3v_f(ji,jj,jk)
[12340]225               END_3D
[6140]226               !
[2528]227            ENDIF
228            !
[3]229         ENDIF
[2528]230         !
[6140]231         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
[11050]232            ! Revert filtered "now" velocities to time split estimate
[12581]233            ! Doing it here also means that asselin filter contribution is removed
[11050]234            zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
[12581]235            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
[4990]236            DO jk = 2, jpkm1
[11050]237               zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
[12581]238               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
[4370]239            END DO
240            DO jk = 1, jpkm1
[11050]241               puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) - (zue(:,:) * r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)) * umask(:,:,jk)
242               pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) - (zve(:,:) * r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)) * vmask(:,:,jk)
[4292]243            END DO
244         ENDIF
245         !
[11050]246      ENDIF ! neuler /= 0
[4354]247      !
248      ! Set "now" and "before" barotropic velocities for next time step:
249      ! JC: Would be more clever to swap variables than to make a full vertical
250      ! integration
251      !
[6140]252      IF(.NOT.ln_linssh ) THEN
[11050]253         hu(:,:,Kmm) = pe3u(:,:,1,Kmm ) * umask(:,:,1)
254         hv(:,:,Kmm) = pe3v(:,:,1,Kmm ) * vmask(:,:,1)
[6140]255         DO jk = 2, jpkm1
[11050]256            hu(:,:,Kmm) = hu(:,:,Kmm) + pe3u(:,:,jk,Kmm ) * umask(:,:,jk)
257            hv(:,:,Kmm) = hv(:,:,Kmm) + pe3v(:,:,jk,Kmm ) * vmask(:,:,jk)
[4354]258         END DO
[11050]259         r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
260         r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
[4354]261      ENDIF
262      !
[11050]263      uu_b(:,:,Kaa) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
264      uu_b(:,:,Kmm) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
265      vv_b(:,:,Kaa) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
266      vv_b(:,:,Kmm) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
[6140]267      DO jk = 2, jpkm1
[11050]268         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
269         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
270         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
271         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
[4354]272      END DO
[11050]273      uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa)
274      vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa)
275      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
276      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
[4354]277      !
[6140]278      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
[11050]279         CALL iom_put(  "ubar", uu_b(:,:,Kmm) )
280         CALL iom_put(  "vbar", vv_b(:,:,Kmm) )
[6140]281      ENDIF
[4990]282      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
[11050]283         zua(:,:,:) = ( puu(:,:,:,Kmm) - zua(:,:,:) ) * z1_2dt
284         zva(:,:,:) = ( pvv(:,:,:,Kmm) - zva(:,:,:) ) * z1_2dt
[10946]285         CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm )
[4990]286      ENDIF
287      !
[12236]288      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   &
289         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask )
[12581]290      !
[9019]291      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
292      IF( l_trddyn     )   DEALLOCATE( zua, zva )
[12581]293      IF( ln_timing    )   CALL timing_stop('dyn_atf_lf')
[2715]294      !
[12581]295   END SUBROUTINE dyn_atf_lf
[3]296
[1502]297   !!=========================================================================
[12581]298END MODULE dynatfLF
Note: See TracBrowser for help on using the repository browser.