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_r14273_HPC-02_Daley_Tiling/src/OCE/DYN – NEMO

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

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

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