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/trunk/src/OCE/DYN – NEMO

source: NEMO/trunk/src/OCE/DYN/dynadv_ubs.F90 @ 13286

Last change on this file since 13286 was 13237, checked in by smasson, 4 years ago

trunk: Mid-year merge, merge back KERNEL-06_techene_e3

  • Property svn:keywords set to Id
File size: 12.7 KB
RevLine 
[643]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   !!======================================================================
[1566]7   !! History :  2.0  ! 2006-08  (R. Benshila, L. Debreu)  Original code
8   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
[643]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
[5062]18   USE trd_oce        ! trends: ocean variables
19   USE trddyn         ! trend manager: dynamics
20   !
[2715]21   USE in_out_manager ! I/O manager
[1129]22   USE prtctl         ! Print control
[2715]23   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
24   USE lib_mpp        ! MPP library
[643]25
26   IMPLICIT NONE
27   PRIVATE
28
[3294]29   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS
[4153]30   REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0   2nd order  ; =1/32 4th order centred
[643]31
[1566]32   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
[643]33
34   !! * Substitutions
[12377]35#  include "do_loop_substitute.h90"
[13237]36#  include "domzgr_substitute.h90"
[643]37   !!----------------------------------------------------------------------
[9598]38   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]39   !! $Id$
[10068]40   !! Software governed by the CeCILL license (see ./LICENSE)
[643]41   !!----------------------------------------------------------------------
42CONTAINS
43
[12377]44   SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs )
[643]45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE dyn_adv_ubs  ***
47      !!
48      !! ** Purpose :   Compute the now momentum advection trend in flux form
[1566]49      !!              and the general trend of the momentum equation.
[643]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
[4153]58      !!                       = 1/32 4th order finite differencing
[643]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
[4153]65      !!      gamma1=1/3 and gamma2=1/32.
[643]66      !!
[12377]67      !! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the 3D advective momentum trends
[643]68      !!
69      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
70      !!----------------------------------------------------------------------
[12377]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
[2715]74      !
[6140]75      INTEGER  ::   ji, jj, jk   ! dummy loop indices
76      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! local scalars
[9019]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
[643]81      !!----------------------------------------------------------------------
[3294]82      !
[643]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
[3294]88      !
[2715]89      zfu_t(:,:,:) = 0._wp
90      zfv_t(:,:,:) = 0._wp
91      zfu_f(:,:,:) = 0._wp
92      zfv_f(:,:,:) = 0._wp
[1566]93      !
[2715]94      zlu_uu(:,:,:,:) = 0._wp
95      zlv_vv(:,:,:,:) = 0._wp 
96      zlu_uv(:,:,:,:) = 0._wp 
97      zlv_vu(:,:,:,:) = 0._wp 
[6140]98      !
99      IF( l_trddyn ) THEN           ! trends: store the input trends
[12377]100         zfu_uw(:,:,:) = puu(:,:,:,Krhs)
101         zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
[1129]102      ENDIF
[1566]103      !                                      ! =========================== !
104      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
105         !                                   ! =========================== !
106         !                                         ! horizontal volume fluxes
[12377]107         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
108         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
[1566]109         !           
[12377]110         DO_2D_00_00
111            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)
112            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)
113            zlu_uv(ji,jj,jk,1) = ( puu (ji  ,jj+1,jk,Kbb) - puu (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
114               &               - ( puu (ji  ,jj  ,jk,Kbb) - puu (ji  ,jj-1,jk,Kbb) ) * fmask(ji  ,jj-1,jk)
115            zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj  ,jk,Kbb) - pvv (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
116               &               - ( pvv (ji  ,jj  ,jk,Kbb) - pvv (ji-1,jj  ,jk,Kbb) ) * fmask(ji-1,jj  ,jk)
117            !
118            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)
119            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)
120            zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
121               &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
122            zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
123               &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
124         END_2D
[1566]125      END DO
[13226]126      CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1.0_wp , zlu_uv(:,:,:,1), 'U', 1.0_wp,  &
127                      &   zlu_uu(:,:,:,2), 'U', 1.0_wp , zlu_uv(:,:,:,2), 'U', 1.0_wp,  & 
128                      &   zlv_vv(:,:,:,1), 'V', 1.0_wp , zlv_vu(:,:,:,1), 'V', 1.0_wp,  &
129                      &   zlv_vv(:,:,:,2), 'V', 1.0_wp , zlv_vu(:,:,:,2), 'V', 1.0_wp   )
[6140]130      !
[1566]131      !                                      ! ====================== !
132      !                                      !  Horizontal advection  !
133      DO jk = 1, jpkm1                       ! ====================== !
134         !                                         ! horizontal volume fluxes
[12377]135         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
136         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
[1566]137         !
[12377]138         DO_2D_10_10
139            zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
140            zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
141            !
142            IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
143            ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
144            ENDIF
145            IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
146            ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
147            ENDIF
148            !
149            zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
150               &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
151               &                * ( zui - gamma1 * zl_u)
152            zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
153               &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
154               &                * ( zvj - gamma1 * zl_v)
155            !
156            zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
157            zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
158            IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
159            ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
160            ENDIF
161            IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
162            ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
163            ENDIF
164            !
165            zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
166               &                * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) - gamma1 * zl_u )
167            zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
168               &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
169         END_2D
170         DO_2D_00_00
171            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
[13237]172               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
173               &                           / e3u(ji,jj,jk,Kmm)
[12377]174            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
[13237]175               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj)   &
176               &                           / e3v(ji,jj,jk,Kmm)
[12377]177         END_2D
[1566]178      END DO
[6140]179      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
[12377]180         zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
181         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
182         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
183         zfu_t(:,:,:) = puu(:,:,:,Krhs)
184         zfv_t(:,:,:) = pvv(:,:,:,Krhs)
[1129]185      ENDIF
[1566]186      !                                      ! ==================== !
187      !                                      !  Vertical advection  !
[6140]188      !                                      ! ==================== !
[12377]189      DO_2D_00_00
190         zfu_uw(ji,jj,jpk) = 0._wp
191         zfv_vw(ji,jj,jpk) = 0._wp
192         zfu_uw(ji,jj, 1 ) = 0._wp
193         zfv_vw(ji,jj, 1 ) = 0._wp
194      END_2D
[6140]195      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
[12377]196         DO_2D_00_00
197            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)
198            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)
199         END_2D
[6140]200      ENDIF
201      DO jk = 2, jpkm1                          ! interior fluxes
[12377]202         DO_2D_01_01
203            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
204         END_2D
205         DO_2D_00_00
206            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) )
207            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) )
208         END_2D
[643]209      END DO
[12377]210      DO_3D_00_00( 1, jpkm1 )
[13237]211         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)   &
212            &                                       / e3u(ji,jj,jk,Kmm)
213         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)   &
214            &                                       / e3v(ji,jj,jk,Kmm)
[12377]215      END_3D
[1566]216      !
[6140]217      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
[12377]218         zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
219         zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
220         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
[1129]221      ENDIF
[6140]222      !                                         ! Control print
[12377]223      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
224         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]225      !
[643]226   END SUBROUTINE dyn_adv_ubs
227
228   !!==============================================================================
229END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.