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 @ 14618

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

#2506 comments and loop optimisation from gm
  • Property svn:keywords set to Id
File size: 15.2 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(jpi,jpj,jpk)   ::   zfu_t, zfu_f, zfu_uw, zfu
83      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f, zfv_vw, zfv, zfw
84      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv
85      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu
86      REAL(wp), DIMENSION(:,:,:), POINTER ::   zpt_u, zpt_v, zpt_w
87      !!----------------------------------------------------------------------
88      !
89      IF( kt == nit000 ) THEN
90         IF(lwp) WRITE(numout,*)
91         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
92         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
93      ENDIF
94      !
95      zfu_t(:,:,:) = 0._wp
96      zfv_t(:,:,:) = 0._wp
97      zfu_f(:,:,:) = 0._wp
98      zfv_f(:,:,:) = 0._wp
99      !
100      zlu_uu(:,:,:,:) = 0._wp
101      zlv_vv(:,:,:,:) = 0._wp 
102      zlu_uv(:,:,:,:) = 0._wp 
103      zlv_vu(:,:,:,:) = 0._wp 
104      !
105      IF( l_trddyn ) THEN           ! trends: store the input trends
106         zfu_uw(:,:,:) = puu(:,:,:,Krhs)
107         zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
108      ENDIF
109      !
110      IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
111         zpt_u => pau(:,:,:)
112         zpt_v => pav(:,:,:)
113         zpt_w => paw(:,:,:)
114      ELSE                          ! MLF: advective velocity = (puu,pvv,ww)
115         zpt_u => puu(:,:,:,Kmm)
116         zpt_v => pvv(:,:,:,Kmm)
117         zpt_w => ww (:,:,:    )
118      ENDIF
119      !
120      !                                      ! =========================== !
121      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
122         !                                   ! =========================== !
123         !                                         ! horizontal volume fluxes
124         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * zpt_u(:,:,jk)
125         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * zpt_v(:,:,jk)
126         !           
127         DO_2D( 0, 0, 0, 0 )                       ! laplacian
128            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)
129            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)
130            zlu_uv(ji,jj,jk,1) = ( puu (ji  ,jj+1,jk,Kbb) - puu (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
131               &               - ( puu (ji  ,jj  ,jk,Kbb) - puu (ji  ,jj-1,jk,Kbb) ) * fmask(ji  ,jj-1,jk)
132            zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj  ,jk,Kbb) - pvv (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
133               &               - ( pvv (ji  ,jj  ,jk,Kbb) - pvv (ji-1,jj  ,jk,Kbb) ) * fmask(ji-1,jj  ,jk)
134            !
135            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)
136            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)
137            zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
138               &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
139            zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
140               &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
141         END_2D
142      END DO
143      CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1.0_wp , zlu_uv(:,:,:,1), 'U', 1.0_wp,  &
144                      &   zlu_uu(:,:,:,2), 'U', 1.0_wp , zlu_uv(:,:,:,2), 'U', 1.0_wp,  & 
145                      &   zlv_vv(:,:,:,1), 'V', 1.0_wp , zlv_vu(:,:,:,1), 'V', 1.0_wp,  &
146                      &   zlv_vv(:,:,:,2), 'V', 1.0_wp , zlv_vu(:,:,:,2), 'V', 1.0_wp   )
147      !
148      !                                      ! ====================== !
149      !                                      !  Horizontal advection  !
150      DO jk = 1, jpkm1                       ! ====================== !
151         !                                         ! horizontal volume fluxes
152         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * zpt_u(:,:,jk)
153         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * zpt_v(:,:,jk)
154         !
155         DO_2D( 1, 0, 1, 0 )                       ! horizontal momentum fluxes at T- and F-point
156            zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
157            zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
158            !
159            IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
160            ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
161            ENDIF
162            IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
163            ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
164            ENDIF
165            !
166            zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
167               &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
168               &                * ( zui - gamma1 * zl_u)
169            zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
170               &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
171               &                * ( zvj - gamma1 * zl_v)
172            !
173            zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
174            zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
175            IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
176            ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
177            ENDIF
178            IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
179            ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
180            ENDIF
181            !
182            zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
183               &                * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) - gamma1 * zl_u )
184            zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
185               &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
186         END_2D
187         DO_2D( 0, 0, 0, 0 )                       ! divergence of horizontal momentum fluxes
188            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
189               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
190               &                           / e3u(ji,jj,jk,Kmm)
191            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
192               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj)   &
193               &                           / e3v(ji,jj,jk,Kmm)
194         END_2D
195      END DO
196      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
197         zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
198         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
199         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
200         zfu_t(:,:,:) = puu(:,:,:,Krhs)
201         zfv_t(:,:,:) = pvv(:,:,:,Krhs)
202      ENDIF
203      !                                      ! ==================== !
204      !                                      !  Vertical advection  !
205      !                                      ! ==================== !
206      !
207      !                                      ! ======================== !
208      IF( PRESENT( no_zad ) ) THEN           !  No vertical advection   !   (except if linear free surface)
209         !                                   ! ======================== !    ------
210         !
211         IF( ln_linssh ) THEN                      ! linear free surface: advection through the surface z=0
212            DO_2D( 0, 0, 0, 0 )
213               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)
214               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)
215               puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj)   &
216                  &                                        / e3u(ji,jj,1,Kmm)
217               pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj)   &
218                  &                                        / e3v(ji,jj,1,Kmm)
219            END_2D
220         ENDIF
221         !                                   ! =================== !
222      ELSE                                   !  Vertical advection !
223         !                                   ! =================== !
224         DO_2D( 0, 0, 0, 0 )                          ! surface/bottom advective fluxes set to zero
225            zfu_uw(ji,jj,jpk) = 0._wp
226            zfv_vw(ji,jj,jpk) = 0._wp
227            zfu_uw(ji,jj, 1 ) = 0._wp
228            zfv_vw(ji,jj, 1 ) = 0._wp
229         END_2D
230         IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
231            DO_2D( 0, 0, 0, 0 )
232               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)
233               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)
234            END_2D
235         ENDIF
236         DO jk = 2, jpkm1                          ! interior fluxes
237            DO_2D( 0, 1, 0, 1 )
238               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
239            END_2D
240            DO_2D( 0, 0, 0, 0 )
241               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) )
242               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) )
243            END_2D
244         END DO
245         DO_3D( 0, 0, 0, 0, 1, jpkm1 )             ! divergence of vertical momentum flux divergence
246            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)   &
247               &                                       / e3u(ji,jj,jk,Kmm)
248            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)   &
249               &                                       / e3v(ji,jj,jk,Kmm)
250         END_3D
251         !
252         IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
253            zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
254            zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
255            CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
256         ENDIF
257         !                                         ! Control print
258         IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
259            &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
260         !
261      ENDIF
262      !
263   END SUBROUTINE dyn_adv_ubs
264
265   !!==============================================================================
266END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.