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.
dynatf_qco.F90 in NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN – NEMO

source: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/DYN/dynatf_qco.F90 @ 14576

Last change on this file since 14576 was 14574, checked in by hadcv, 3 years ago

#2600: Merge in trunk changes to r14509
  • Property svn:keywords set to Id
File size: 15.2 KB
Line 
1MODULE dynatf_qco
2   !!=========================================================================
3   !!                       ***  MODULE  dynatf_qco  ***
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 -> dynatfLF.F90. Now just does time filtering.
23   !!-------------------------------------------------------------------------
24
25   !!----------------------------------------------------------------------------------------------
26   !!   dyn_atf_qco       : 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   USE zdfdrg    , ONLY: ln_drgice_imp, rCdU_top
46   !
47   USE in_out_manager ! I/O manager
48   USE iom            ! I/O manager library
49   USE lbclnk         ! lateral boundary condition (or mpp link)
50   USE lib_mpp        ! MPP library
51   USE prtctl         ! Print control
52   USE timing         ! Timing
53
54   IMPLICIT NONE
55   PRIVATE
56
57   PUBLIC    dyn_atf_qco   ! routine called by step.F90
58
59   !! * Substitutions
60#  include "do_loop_substitute.h90"
61#  include "domzgr_substitute.h90"
62   !!----------------------------------------------------------------------
63   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
64   !! $Id$
65   !! Software governed by the CeCILL license (see ./LICENSE)
66   !!----------------------------------------------------------------------
67CONTAINS
68
69   SUBROUTINE dyn_atf_qco( kt, Kbb, Kmm, Kaa, puu, pvv )
70      !!----------------------------------------------------------------------
71      !!                  ***  ROUTINE dyn_atf_qco  ***
72      !!
73      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
74      !!             condition on the after velocity and apply the Asselin time
75      !!             filter to the now fields.
76      !!
77      !! ** Method  : * Ensure after velocities transport matches time splitting
78      !!             estimate (ln_dynspg_ts=T)
79      !!
80      !!              * Apply lateral boundary conditions on after velocity
81      !!             at the local domain boundaries through lbc_lnk call,
82      !!             at the one-way open boundaries (ln_bdy=T),
83      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
84      !!
85      !!              * Apply the Asselin time filter to the now fields
86      !!             arrays to start the next time step:
87      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm))
88      !!                                    + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ]
89      !!             Note that with flux form advection and non linear free surface,
90      !!             the time filter is applied on thickness weighted velocity.
91      !!             As a result, dyn_atf_lf MUST be called after tra_atf.
92      !!
93      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity
94      !!----------------------------------------------------------------------
95      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
96      INTEGER                             , INTENT(in   ) :: Kbb, Kmm, Kaa    ! before and after time level indices
97      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered
98      !
99      INTEGER  ::   ji, jj, jk   ! dummy loop indices
100      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars
101      REAL(wp) ::   zve3a, zve3n, zve3b           !   -      -
102      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
103      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zua, zva
104      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zutau, zvtau
105      !!----------------------------------------------------------------------
106      !
107      IF( ln_timing    )   CALL timing_start('dyn_atf_qco')
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_atf_qco : Asselin time filtering'
114         IF(lwp) WRITE(numout,*) '~~~~~~~'
115      ENDIF
116      !
117      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
118         !
119         !                                  ! Kinetic energy and Conversion
120         IF( ln_KE_trd  )   CALL trd_dyn( puu(:,:,:,Kaa), pvv(:,:,:,Kaa), jpdyn_ken, kt, Kmm )
121         !
122         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
123            zua(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) * r1_Dt
124            zva(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) * r1_Dt
125            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
126            CALL iom_put( "vtrd_tot", zva )
127         ENDIF
128         !
129         zua(:,:,:) = puu(:,:,:,Kmm)             ! save the now velocity before the asselin filter
130         zva(:,:,:) = pvv(:,:,:,Kmm)             ! (caution: there will be a shift by 1 timestep in the
131         !                                  !  computation of the asselin filter trends)
132      ENDIF
133
134      ! Time filter and swap of dynamics arrays
135      ! ------------------------------------------
136
137      IF( .NOT. l_1st_euler ) THEN    !* Leap-Frog : Asselin time filter
138         !                                ! =============!
139         IF( ln_linssh ) THEN             ! Fixed volume !
140            !                             ! =============!
141            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
142               puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
143               pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
144            END_3D
145            !                             ! ================!
146         ELSE                             ! Variable volume !
147            !                             ! ================!
148            !
149            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
150               ! Before filtered scale factor at (u/v)-points
151               DO_3D( 1, 1, 1, 1, 1, jpkm1 )
152                  puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
153                  pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
154               END_3D
155               !
156            ELSE                          ! Asselin filter applied on thickness weighted velocity
157               !
158               DO_3D( 1, 1, 1, 1, 1, jpkm1 )
159                  zue3a = ( 1._wp + r3u(ji,jj,Kaa) * umask(ji,jj,jk) ) * puu(ji,jj,jk,Kaa)
160                  zve3a = ( 1._wp + r3v(ji,jj,Kaa) * vmask(ji,jj,jk) ) * pvv(ji,jj,jk,Kaa)
161                  zue3n = ( 1._wp + r3u(ji,jj,Kmm) * umask(ji,jj,jk) ) * puu(ji,jj,jk,Kmm)
162                  zve3n = ( 1._wp + r3v(ji,jj,Kmm) * vmask(ji,jj,jk) ) * pvv(ji,jj,jk,Kmm)
163                  zue3b = ( 1._wp + r3u(ji,jj,Kbb) * umask(ji,jj,jk) ) * puu(ji,jj,jk,Kbb)
164                  zve3b = ( 1._wp + r3v(ji,jj,Kbb) * vmask(ji,jj,jk) ) * pvv(ji,jj,jk,Kbb)
165                  !                                                 ! filtered scale factor at U-,V-points
166                  puu(ji,jj,jk,Kmm) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ( 1._wp + r3u_f(ji,jj)*umask(ji,jj,jk) )
167                  pvv(ji,jj,jk,Kmm) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ( 1._wp + r3v_f(ji,jj)*vmask(ji,jj,jk) )
168               END_3D
169               !
170            ENDIF
171            !
172         ENDIF
173         !
174         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
175            ! Revert filtered "now" velocities to time split estimate
176            ! Doing it here also means that asselin filter contribution is removed
177            ! zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
178            ! zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
179            ! DO jk = 2, jpkm1
180            !    zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
181            !    zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
182            ! END DO
183            zue(:,:) = e3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
184            zve(:,:) = e3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
185            DO jk = 2, jpkm1
186               zue(:,:) = zue(:,:) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
187               zve(:,:) = zve(:,:) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
188            END DO
189            DO jk = 1, jpkm1
190               puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) - (zue(:,:) * r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)) * umask(:,:,jk)
191               pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) - (zve(:,:) * r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)) * vmask(:,:,jk)
192            END DO
193         ENDIF
194         !
195      ENDIF ! .NOT. l_1st_euler
196      !
197      ! Set "now" and "before" barotropic velocities for next time step:
198      ! JC: Would be more clever to swap variables than to make a full vertical
199      ! integration
200      ! CAUTION : calculation need to be done in the same way than see GM
201#if defined key_linssh
202      uu_b(:,:,Kaa) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
203      uu_b(:,:,Kmm) = e3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
204      vv_b(:,:,Kaa) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
205      vv_b(:,:,Kmm) = e3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
206      DO jk = 2, jpkm1
207         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
208         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
209         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
210         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
211      END DO
212      uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa)
213      vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa)
214      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
215      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
216#else
217      uu_b(:,:,Kaa) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
218      uu_b(:,:,Kmm) = (e3u_0(:,:,1) * ( 1._wp + r3u_f(:,:) * umask(:,:,1) )) * puu(:,:,1,Kmm) * umask(:,:,1)
219      vv_b(:,:,Kaa) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
220      vv_b(:,:,Kmm) = (e3v_0(:,:,1) * ( 1._wp + r3v_f(:,:) * vmask(:,:,1))) * pvv(:,:,1,Kmm) * vmask(:,:,1)
221      DO jk = 2, jpkm1
222         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
223         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + (e3u_0(:,:,jk) * ( 1._wp + r3u_f(:,:) * umask(:,:,jk) )) * puu(:,:,jk,Kmm) * umask(:,:,jk)
224         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
225         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + (e3v_0(:,:,jk) * ( 1._wp + r3v_f(:,:) * vmask(:,:,jk) )) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
226      END DO
227      uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa)
228      vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa)
229      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * (r1_hu_0(:,:)/( 1._wp + r3u_f(:,:) ))
230      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * (r1_hv_0(:,:)/( 1._wp + r3v_f(:,:) ))
231#endif
232      !
233      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
234         CALL iom_put(  "ubar", uu_b(:,:,Kmm) )
235         CALL iom_put(  "vbar", vv_b(:,:,Kmm) )
236      ENDIF
237      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
238         zua(:,:,:) = ( puu(:,:,:,Kmm) - zua(:,:,:) ) * r1_Dt
239         zva(:,:,:) = ( pvv(:,:,:,Kmm) - zva(:,:,:) ) * r1_Dt
240         CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm )
241      ENDIF
242      !
243      IF ( iom_use("utau") ) THEN
244         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
245            ALLOCATE(zutau(jpi,jpj))
246            DO_2D( 0, 0, 0, 0 )
247               jk = miku(ji,jj)
248               zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
249            END_2D
250            CALL iom_put(  "utau", zutau(:,:) )
251            DEALLOCATE(zutau)
252         ELSE
253            CALL iom_put(  "utau", utau(:,:) )
254         ENDIF
255      ENDIF
256      !
257      IF ( iom_use("vtau") ) THEN
258         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
259            ALLOCATE(zvtau(jpi,jpj))
260            DO_2D( 0, 0, 0, 0 )
261               jk = mikv(ji,jj)
262               zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
263            END_2D
264            CALL iom_put(  "vtau", zvtau(:,:) )
265            DEALLOCATE(zvtau)
266         ELSE
267            CALL iom_put(  "vtau", vtau(:,:) )
268         ENDIF
269      ENDIF
270      !
271      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   &
272         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask )
273      !
274      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
275      IF( l_trddyn     )   DEALLOCATE( zua, zva )
276      IF( ln_timing    )   CALL timing_stop('dyn_atf_qco')
277      !
278   END SUBROUTINE dyn_atf_qco
279
280   !!=========================================================================
281END MODULE dynatf_qco
Note: See TracBrowser for help on using the repository browser.