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.
dynadv_ubs.F90 in NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN – NEMO

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynadv_ubs.F90 @ 15574

Last change on this file since 15574 was 15574, checked in by techene, 3 years ago

#2605 #2715 trunk merged into dev_r14318_RK3_stage1

  • Property svn:keywords set to Id
File size: 17.4 KB
Line 
1MODULE dynadv_ubs
2   !!======================================================================
3   !!                       ***  MODULE  dynadv_ubs  ***
4   !! Ocean dynamics: Update the momentum trend with the flux form advection
5   !!                 trend using a 3rd order upstream biased scheme
6   !!======================================================================
7   !! History :  2.0  ! 2006-08  (R. Benshila, L. Debreu)  Original code
8   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
9   !!----------------------------------------------------------------------
10
11   !!----------------------------------------------------------------------
12   !!   dyn_adv_ubs   : flux form momentum advection using    (ln_dynadv=T)
13   !!                   an 3rd order Upstream Biased Scheme or Quick scheme
14   !!                   combined with 2nd or 4th order finite differences
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers
17   USE dom_oce        ! ocean space and time domain
18   USE trd_oce        ! trends: ocean variables
19   USE trddyn         ! trend manager: dynamics
20   !
21   USE in_out_manager ! I/O manager
22   USE prtctl         ! Print control
23   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
24   USE lib_mpp        ! MPP library
25
26   IMPLICIT NONE
27   PRIVATE
28
29   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS
30   REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0   2nd order  ; =1/32 4th order centred
31
32   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
33
34   !! * Substitutions
35#  include "do_loop_substitute.h90"
36#  include "domzgr_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
39   !! $Id$
40   !! Software governed by the CeCILL license (see ./LICENSE)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE dyn_adv_ubs  ***
47      !!
48      !! ** Purpose :   Compute the now momentum advection trend in flux form
49      !!              and the general trend of the momentum equation.
50      !!
51      !! ** Method  :   The scheme is the one implemeted in ROMS. It depends
52      !!      on two parameter gamma1 and gamma2. The former control the
53      !!      upstream baised part of the scheme and the later the centred
54      !!      part:     gamma1 = 0    pure centered  (no diffusive part)
55      !!                       = 1/4  Quick scheme
56      !!                       = 1/3  3rd order Upstream biased scheme
57      !!                gamma2 = 0    2nd order finite differencing
58      !!                       = 1/32 4th order finite differencing
59      !!      For stability reasons, the first term of the fluxes which cor-
60      !!      responds to a second order centered scheme is evaluated using 
61      !!      the now velocity (centered in time) while the second term which 
62      !!      is the diffusive part of the scheme, is evaluated using the
63      !!      before velocity (forward in time).
64      !!      Default value (hard coded in the begining of the module) are
65      !!      gamma1=1/3 and gamma2=1/32.
66      !!
67      !!                In RK3 time stepping case, the optional arguments
68      !!      (pau,pav,paw) are present. They are used as advective velocity 
69      !!      while the advected velocity remains (puu,pvv).
70      !!
71      !! ** Action  :   (puu,pvv)(:,:,:,Krhs)   updated with the advective trend
72      !!
73      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
74      !!----------------------------------------------------------------------
75      INTEGER                                     , INTENT(in   ) ::   kt , Kbb, Kmm, Krhs   ! ocean time-step and level indices
76      INTEGER                   , OPTIONAL        , INTENT(in   ) ::   no_zad                ! no vertical advection compotation
77      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) ::   puu, pvv              ! ocean velocities and RHS of momentum equation
78      REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in   ) ::   pau, pav, paw         ! advective velocity
79      !
80      INTEGER  ::   ji, jj, jk   ! dummy loop indices
81      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v, zzu, zzv   ! local scalars
82      REAL(wp), DIMENSION(A2D(nn_hls),jpk)   ::   zfu_t, zfu_f, zfu_uw, zfu
83      REAL(wp), DIMENSION(A2D(nn_hls),jpk)   ::   zfv_t, zfv_f, zfv_vw, zfv, zfw
84      REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) ::   zlu_uu, zlu_uv
85      REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) ::   zlv_vv, zlv_vu
86      REAL(wp), DIMENSION(:,:,:), POINTER ::   zpt_u, zpt_v, zpt_w
87      !!----------------------------------------------------------------------
88      !
89      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
90         IF( kt == nit000 ) THEN
91            IF(lwp) WRITE(numout,*)
92            IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
93            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
94         ENDIF
95      ENDIF
96      !
97      zfu_t(:,:,:) = 0._wp
98      zfv_t(:,:,:) = 0._wp
99      zfu_f(:,:,:) = 0._wp
100      zfv_f(:,:,:) = 0._wp
101      !
102      zlu_uu(:,:,:,:) = 0._wp
103      zlv_vv(:,:,:,:) = 0._wp 
104      zlu_uv(:,:,:,:) = 0._wp 
105      zlv_vu(:,:,:,:) = 0._wp 
106      !
107      IF( l_trddyn ) THEN           ! trends: store the input trends
108         zfu_uw(:,:,:) = puu(:,:,:,Krhs)
109         zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
110      ENDIF
111      !
112      IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
113         zpt_u => pau(:,:,:)
114         zpt_v => pav(:,:,:)
115         zpt_w => paw(:,:,:)
116      ELSE                          ! MLF: advective velocity = (puu,pvv,ww)
117         zpt_u => puu(:,:,:,Kmm)
118         zpt_v => pvv(:,:,:,Kmm)
119         zpt_w => ww (:,:,:    )
120      ENDIF
121      !
122      !                                      ! =========================== !
123      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
124         !                                   ! =========================== !
125         !                                         ! horizontal volume fluxes
126         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
127            zfu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
128            zfv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
129         END_2D
130         !           
131         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                       ! laplacian
132            ! round brackets added to fix the order of floating point operations
133            ! needed to ensure halo 1 - halo 2 compatibility
134            zlu_uu(ji,jj,jk,1) = ( ( puu (ji+1,jj  ,jk,Kbb) - puu (ji  ,jj  ,jk,Kbb) &
135               &                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
136               &                 + ( puu (ji-1,jj  ,jk,Kbb) - puu (ji  ,jj  ,jk,Kbb) &
137               &                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
138               &                 ) * umask(ji  ,jj  ,jk)
139            zlv_vv(ji,jj,jk,1) = ( ( pvv (ji  ,jj+1,jk,Kbb) - pvv (ji  ,jj  ,jk,Kbb) &
140               &                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
141               &                 + ( pvv (ji  ,jj-1,jk,Kbb) - pvv (ji  ,jj  ,jk,Kbb) &
142               &                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
143               &                 ) * vmask(ji  ,jj  ,jk)
144            zlu_uv(ji,jj,jk,1) = (  puu (ji  ,jj+1,jk,Kbb) - puu (ji  ,jj  ,jk,Kbb)  ) * fmask(ji  ,jj  ,jk)   &
145               &               - (  puu (ji  ,jj  ,jk,Kbb) - puu (ji  ,jj-1,jk,Kbb)  ) * fmask(ji  ,jj-1,jk)
146            zlv_vu(ji,jj,jk,1) = (  pvv (ji+1,jj  ,jk,Kbb) - pvv (ji  ,jj  ,jk,Kbb)  ) * fmask(ji  ,jj  ,jk)   &
147               &               - (  pvv (ji  ,jj  ,jk,Kbb) - pvv (ji-1,jj  ,jk,Kbb)  ) * fmask(ji-1,jj  ,jk)
148            !
149            ! round brackets added to fix the order of floating point operations
150            ! needed to ensure halo 1 - halo 2 compatibility
151            zlu_uu(ji,jj,jk,2) = ( ( zfu(ji+1,jj  ,jk) - zfu(ji  ,jj  ,jk)           &
152               &                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
153               &                 + ( zfu(ji-1,jj  ,jk) - zfu(ji  ,jj  ,jk)           &
154               &                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
155               &                 ) * umask(ji  ,jj  ,jk)
156            zlv_vv(ji,jj,jk,2) = ( ( zfv(ji  ,jj+1,jk) - zfv(ji  ,jj  ,jk)           &
157               &                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
158               &                 + ( zfv(ji  ,jj-1,jk) - zfv(ji  ,jj  ,jk)           &
159               &                   )                                                 & ! bracket for halo 1 - halo 2 compatibility
160               &                 ) * vmask(ji  ,jj  ,jk)
161            zlu_uv(ji,jj,jk,2) = (  zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk)  ) * fmask(ji  ,jj  ,jk)             &
162               &               - (  zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk)  ) * fmask(ji  ,jj-1,jk)
163            zlv_vu(ji,jj,jk,2) = (  zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk)  ) * fmask(ji  ,jj  ,jk)             &
164               &               - (  zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk)  ) * fmask(ji-1,jj  ,jk)
165         END_2D
166      END DO
167      IF( nn_hls == 1 ) CALL lbc_lnk( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', -1.0_wp , zlu_uv(:,:,:,1), 'U', -1.0_wp,  &
168                                              &   zlu_uu(:,:,:,2), 'U', -1.0_wp , zlu_uv(:,:,:,2), 'U', -1.0_wp,  &
169                                              &   zlv_vv(:,:,:,1), 'V', -1.0_wp , zlv_vu(:,:,:,1), 'V', -1.0_wp,  &
170                                              &   zlv_vv(:,:,:,2), 'V', -1.0_wp , zlv_vu(:,:,:,2), 'V', -1.0_wp   )
171      !
172      !                                      ! ====================== !
173      !                                      !  Horizontal advection  !
174      DO jk = 1, jpkm1                       ! ====================== !
175         !                                         ! horizontal volume fluxes
176         DO_2D( 1, 1, 1, 1 )
177            zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
178            zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
179         END_2D
180         !
181         DO_2D( 1, 0, 1, 0 )                       ! horizontal momentum fluxes at T- and F-point
182            zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
183            zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
184            !
185            IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
186            ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
187            ENDIF
188            IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
189            ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
190            ENDIF
191            !
192            zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
193               &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
194               &                * ( zui - gamma1 * zl_u)
195            zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
196               &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
197               &                * ( zvj - gamma1 * zl_v)
198            !
199            zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
200            zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
201            IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
202            ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
203            ENDIF
204            IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
205            ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
206            ENDIF
207            !
208            zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
209               &                * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) - gamma1 * zl_u )
210            zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
211               &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
212         END_2D
213         DO_2D( 0, 0, 0, 0 )                       ! divergence of horizontal momentum fluxes
214            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
215               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
216               &                           / e3u(ji,jj,jk,Kmm)
217            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
218               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj)   &
219               &                           / e3v(ji,jj,jk,Kmm)
220         END_2D
221      END DO
222      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
223         zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
224         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
225         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
226         zfu_t(:,:,:) = puu(:,:,:,Krhs)
227         zfv_t(:,:,:) = pvv(:,:,:,Krhs)
228      ENDIF
229      !                                      ! ==================== !
230      !                                      !  Vertical advection  !
231      !                                      ! ==================== !
232      !
233      !                                      ! ======================== !
234      IF( PRESENT( no_zad ) ) THEN           !  No vertical advection   !   (except if linear free surface)
235         !                                   ! ======================== !    ------
236         !
237         IF( ln_linssh ) THEN                      ! linear free surface: advection through the surface z=0
238            DO_2D( 0, 0, 0, 0 )
239               zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
240               zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
241               puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj)   &
242                  &                                        / e3u(ji,jj,1,Kmm)
243               pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj)   &
244                  &                                        / e3v(ji,jj,1,Kmm)
245            END_2D
246         ENDIF
247         !                                   ! =================== !
248      ELSE                                   !  Vertical advection !
249         !                                   ! =================== !
250         DO_2D( 0, 0, 0, 0 )                          ! surface/bottom advective fluxes set to zero
251            zfu_uw(ji,jj,jpk) = 0._wp
252            zfv_vw(ji,jj,jpk) = 0._wp
253            zfu_uw(ji,jj, 1 ) = 0._wp
254            zfv_vw(ji,jj, 1 ) = 0._wp
255         END_2D
256         IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
257            DO_2D( 0, 0, 0, 0 )
258               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
259               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
260            END_2D
261         ENDIF
262         DO jk = 2, jpkm1                          ! interior fluxes
263            DO_2D( 0, 1, 0, 1 )
264               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
265            END_2D
266            DO_2D( 0, 0, 0, 0 )
267               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
268               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
269            END_2D
270         END DO
271         DO_3D( 0, 0, 0, 0, 1, jpkm1 )             ! divergence of vertical momentum flux divergence
272            puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj)   &
273               &                                       / e3u(ji,jj,jk,Kmm)
274            pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj)   &
275               &                                       / e3v(ji,jj,jk,Kmm)
276         END_3D
277         !
278         IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
279            zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
280            zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
281            CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
282         ENDIF
283         !                                         ! Control print
284         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
285            &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
286         !
287      ENDIF
288      !
289   END SUBROUTINE dyn_adv_ubs
290
291   !!==============================================================================
292END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.