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

Last change on this file since 14381 was 14381, checked in by techene, 8 months ago

#2605 : add advection velocity as input in advective flux form routines

  • Property svn:keywords set to Id
File size: 13.7 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 )
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      !! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the 3D advective momentum trends
68      !!
69      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
70      !!----------------------------------------------------------------------
71      INTEGER                                               , INTENT(in   ) ::   kt              ! ocean time-step index
72      INTEGER                                               , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
73      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt)          , TARGET, INTENT(inout) ::   puu, pvv        ! ocean velocities and RHS of momentum equation
74      REAL(wp), DIMENSION(:,:,:)          , OPTIONAL, TARGET, INTENT(in   ) ::   pau, pav, paw   ! advective velocity
75      !
76      INTEGER  ::   ji, jj, jk   ! dummy loop indices
77      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! local scalars
78      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f, zfu_uw, zfu
79      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f, zfv_vw, zfv, zfw
80      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv
81      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu
82      REAL(wp), DIMENSION(:,:,:), POINTER ::   zptu, zptv, zptw
83      !!----------------------------------------------------------------------
84      !
85      IF( kt == nit000 ) THEN
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
88         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
89      ENDIF
90      !
91      zfu_t(:,:,:) = 0._wp
92      zfv_t(:,:,:) = 0._wp
93      zfu_f(:,:,:) = 0._wp
94      zfv_f(:,:,:) = 0._wp
95      !
96      zlu_uu(:,:,:,:) = 0._wp
97      zlv_vv(:,:,:,:) = 0._wp 
98      zlu_uv(:,:,:,:) = 0._wp 
99      zlv_vu(:,:,:,:) = 0._wp 
100      !
101      IF( l_trddyn ) THEN           ! trends: store the input trends
102         zfu_uw(:,:,:) = puu(:,:,:,Krhs)
103         zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
104      ENDIF
105      !
106      IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
107         zptu => pau(:,:,:)
108         zptv => pav(:,:,:)
109         zptw => paw(:,:,:)
110      ELSE                          ! MLF: advective velocity = (puu,pvv,ww)
111         zptu => puu(:,:,:,Kmm)
112         zptv => pvv(:,:,:,Kmm)
113         zptw => ww (:,:,:    )
114      ENDIF
115      !
116      !                                      ! =========================== !
117      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
118         !                                   ! =========================== !
119         !                                         ! horizontal volume fluxes
120         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * zptu(:,:,jk)
121         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * zptv(:,:,jk)
122         !           
123         DO_2D( 0, 0, 0, 0 )                       ! laplacian
124            zlu_uu(ji,jj,jk,1) = ( puu (ji+1,jj  ,jk,Kbb) - 2.*puu (ji,jj,jk,Kbb) + puu (ji-1,jj  ,jk,Kbb) ) * umask(ji,jj,jk)
125            zlv_vv(ji,jj,jk,1) = ( pvv (ji  ,jj+1,jk,Kbb) - 2.*pvv (ji,jj,jk,Kbb) + pvv (ji  ,jj-1,jk,Kbb) ) * vmask(ji,jj,jk)
126            zlu_uv(ji,jj,jk,1) = ( puu (ji  ,jj+1,jk,Kbb) - puu (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
127               &               - ( puu (ji  ,jj  ,jk,Kbb) - puu (ji  ,jj-1,jk,Kbb) ) * fmask(ji  ,jj-1,jk)
128            zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj  ,jk,Kbb) - pvv (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
129               &               - ( pvv (ji  ,jj  ,jk,Kbb) - pvv (ji-1,jj  ,jk,Kbb) ) * fmask(ji-1,jj  ,jk)
130            !
131            zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj  ,jk) - 2.*zfu(ji,jj,jk) + zfu(ji-1,jj  ,jk) ) * umask(ji,jj,jk)
132            zlv_vv(ji,jj,jk,2) = ( zfv(ji  ,jj+1,jk) - 2.*zfv(ji,jj,jk) + zfv(ji  ,jj-1,jk) ) * vmask(ji,jj,jk)
133            zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
134               &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
135            zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
136               &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
137         END_2D
138      END DO
139      CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1.0_wp , zlu_uv(:,:,:,1), 'U', 1.0_wp,  &
140                      &   zlu_uu(:,:,:,2), 'U', 1.0_wp , zlu_uv(:,:,:,2), 'U', 1.0_wp,  & 
141                      &   zlv_vv(:,:,:,1), 'V', 1.0_wp , zlv_vu(:,:,:,1), 'V', 1.0_wp,  &
142                      &   zlv_vv(:,:,:,2), 'V', 1.0_wp , zlv_vu(:,:,:,2), 'V', 1.0_wp   )
143      !
144      !                                      ! ====================== !
145      !                                      !  Horizontal advection  !
146      DO jk = 1, jpkm1                       ! ====================== !
147         !                                         ! horizontal volume fluxes
148         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * zptu(:,:,jk)
149         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * zptv(:,:,jk)
150         !
151         DO_2D( 1, 0, 1, 0 )                       ! horizontal momentum fluxes at T- and F-point
152            zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
153            zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
154            !
155            IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
156            ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
157            ENDIF
158            IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
159            ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
160            ENDIF
161            !
162            zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
163               &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
164               &                * ( zui - gamma1 * zl_u)
165            zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
166               &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
167               &                * ( zvj - gamma1 * zl_v)
168            !
169            zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
170            zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
171            IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
172            ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
173            ENDIF
174            IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
175            ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
176            ENDIF
177            !
178            zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
179               &                * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) - gamma1 * zl_u )
180            zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
181               &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
182         END_2D
183         DO_2D( 0, 0, 0, 0 )                       ! divergence of horizontal momentum fluxes
184            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
185               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
186               &                           / e3u(ji,jj,jk,Kmm)
187            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
188               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj)   &
189               &                           / e3v(ji,jj,jk,Kmm)
190         END_2D
191      END DO
192      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
193         zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
194         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
195         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
196         zfu_t(:,:,:) = puu(:,:,:,Krhs)
197         zfv_t(:,:,:) = pvv(:,:,:,Krhs)
198      ENDIF
199      !                                      ! ==================== !
200      !                                      !  Vertical advection  !
201      !                                      ! ==================== !
202      DO_2D( 0, 0, 0, 0 )                          ! surface/bottom advective fluxes set to zero
203         zfu_uw(ji,jj,jpk) = 0._wp
204         zfv_vw(ji,jj,jpk) = 0._wp
205         zfu_uw(ji,jj, 1 ) = 0._wp
206         zfv_vw(ji,jj, 1 ) = 0._wp
207      END_2D
208      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
209         DO_2D( 0, 0, 0, 0 )
210            zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zptw(ji,jj,1) + e1e2t(ji+1,jj) * zptw(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
211            zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zptw(ji,jj,1) + e1e2t(ji,jj+1) * zptw(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
212         END_2D
213      ENDIF
214      DO jk = 2, jpkm1                          ! interior fluxes
215         DO_2D( 0, 1, 0, 1 )
216            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zptw(ji,jj,jk)
217         END_2D
218         DO_2D( 0, 0, 0, 0 )
219            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) )
220            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) )
221         END_2D
222      END DO
223      DO_3D( 0, 0, 0, 0, 1, jpkm1 )             ! divergence of vertical momentum flux divergence
224         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)   &
225            &                                       / e3u(ji,jj,jk,Kmm)
226         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)   &
227            &                                       / e3v(ji,jj,jk,Kmm)
228      END_3D
229      !
230      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
231         zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
232         zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
233         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
234      ENDIF
235      !                                         ! Control print
236      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
237         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
238      !
239   END SUBROUTINE dyn_adv_ubs
240
241   !!==============================================================================
242END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.