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 branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 8882

Last change on this file since 8882 was 8882, checked in by flavoni, 6 years ago

dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864

  • Property svn:keywords set to Id
File size: 13.1 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   USE timing         ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS
31   REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0   2nd order  ; =1/32 4th order centred
32
33   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
34
35   !! * Substitutions
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
39   !! $Id$
40   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42CONTAINS
43
44   SUBROUTINE dyn_adv_ubs( kt )
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 : - (ua,va) 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      !
73      INTEGER  ::   ji, jj, jk   ! dummy loop indices
74      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! local scalars
75      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f, zfu_uw, zfu
76      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f, zfv_vw, zfv, zfw
77      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv
78      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu
79      !!----------------------------------------------------------------------
80      !
81      IF( ln_timing )   CALL timing_start('dyn_adv_ubs')
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(:,:,:) = ua(:,:,:)
101         zfv_vw(:,:,:) = va(:,:,:)
102      ENDIF
103      !                                      ! =========================== !
104      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
105         !                                   ! =========================== !
106         !                                         ! horizontal volume fluxes
107         zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
108         zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
109         !           
110         DO jj = 2, jpjm1                          ! laplacian
111            DO ji = fs_2, fs_jpim1   ! vector opt.
112               zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj  ,jk) - 2.*ub (ji,jj,jk) + ub (ji-1,jj  ,jk) ) * umask(ji,jj,jk)
113               zlv_vv(ji,jj,jk,1) = ( vb (ji  ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji  ,jj-1,jk) ) * vmask(ji,jj,jk)
114               zlu_uv(ji,jj,jk,1) = ( ub (ji  ,jj+1,jk) - ub (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
115                  &               - ( ub (ji  ,jj  ,jk) - ub (ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
116               zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj  ,jk) - vb (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
117                  &               - ( vb (ji  ,jj  ,jk) - vb (ji-1,jj  ,jk) ) * 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 DO
126         END DO
127      END DO
128      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
129      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
130      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
131      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 
132      !
133      !                                      ! ====================== !
134      !                                      !  Horizontal advection  !
135      DO jk = 1, jpkm1                       ! ====================== !
136         !                                         ! horizontal volume fluxes
137         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
138         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
139         !
140         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
141            DO ji = 1, fs_jpim1   ! vector opt.
142               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
143               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
144               !
145               IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
146               ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
147               ENDIF
148               IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
149               ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
150               ENDIF
151               !
152               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
153                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
154                  &                * ( zui - gamma1 * zl_u)
155               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
156                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
157                  &                * ( zvj - gamma1 * zl_v)
158               !
159               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
160               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
161               IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
162               ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
163               ENDIF
164               IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
165               ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
166               ENDIF
167               !
168               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
169                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
170               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
171                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
172            END DO
173         END DO
174         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
177                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
178               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
179                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
180            END DO
181         END DO
182      END DO
183      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
184         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
185         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
186         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
187         zfu_t(:,:,:) = ua(:,:,:)
188         zfv_t(:,:,:) = va(:,:,:)
189      ENDIF
190      !                                      ! ==================== !
191      !                                      !  Vertical advection  !
192      !                                      ! ==================== !
193      DO jj = 2, jpjm1                             ! surface/bottom advective fluxes set to zero                 
194         DO ji = fs_2, fs_jpim1
195            zfu_uw(ji,jj,jpk) = 0._wp
196            zfv_vw(ji,jj,jpk) = 0._wp
197            zfu_uw(ji,jj, 1 ) = 0._wp
198            zfv_vw(ji,jj, 1 ) = 0._wp
199         END DO
200      END DO
201      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
202         DO jj = 2, jpjm1
203            DO ji = fs_2, fs_jpim1
204               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1)
205               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1)
206            END DO
207         END DO
208      ENDIF
209      DO jk = 2, jpkm1                          ! interior fluxes
210         DO jj = 2, jpj
211            DO ji = 2, jpi
212               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
213            END DO
214         END DO
215         DO jj = 2, jpjm1
216            DO ji = fs_2, fs_jpim1   ! vector opt.
217               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
218               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
219            END DO
220         END DO
221      END DO
222      DO jk = 1, jpkm1                          ! divergence of vertical momentum flux divergence
223         DO jj = 2, jpjm1
224            DO ji = fs_2, fs_jpim1   ! vector opt.
225               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
226               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
227            END DO
228         END DO
229      END DO
230      !
231      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
232         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
233         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
234         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
235      ENDIF
236      !                                         ! Control print
237      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
238         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
239      !
240      IF( ln_timing )   CALL timing_stop('dyn_adv_ubs')
241      !
242   END SUBROUTINE dyn_adv_ubs
243
244   !!==============================================================================
245END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.