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/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 5904

Last change on this file since 5904 was 5883, checked in by gm, 9 years ago

#1613: vvl by default: TRA/TRC remove optimization associated with linear free surface
  • Property svn:keywords set to Id
File size: 13.6 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 wrk_nemo       ! Memory Allocation
26   USE timing         ! Timing
27
28   IMPLICIT NONE
29   PRIVATE
30
31   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS
32   REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0   2nd order  ; =1/32 4th order centred
33
34   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
35
36   !! * Substitutions
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE dyn_adv_ubs( kt )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE dyn_adv_ubs  ***
48      !!
49      !! ** Purpose :   Compute the now momentum advection trend in flux form
50      !!              and the general trend of the momentum equation.
51      !!
52      !! ** Method  :   The scheme is the one implemeted in ROMS. It depends
53      !!      on two parameter gamma1 and gamma2. The former control the
54      !!      upstream baised part of the scheme and the later the centred
55      !!      part:     gamma1 = 0    pure centered  (no diffusive part)
56      !!                       = 1/4  Quick scheme
57      !!                       = 1/3  3rd order Upstream biased scheme
58      !!                gamma2 = 0    2nd order finite differencing
59      !!                       = 1/32 4th order finite differencing
60      !!      For stability reasons, the first term of the fluxes which cor-
61      !!      responds to a second order centered scheme is evaluated using 
62      !!      the now velocity (centered in time) while the second term which 
63      !!      is the diffusive part of the scheme, is evaluated using the
64      !!      before velocity (forward in time).
65      !!      Default value (hard coded in the begining of the module) are
66      !!      gamma1=1/3 and gamma2=1/32.
67      !!
68      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
69      !!
70      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
71      !!----------------------------------------------------------------------
72      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
73      !
74      INTEGER  ::   ji, jj, jk            ! dummy loop indices
75      REAL(wp) ::   zbu, zbv    ! temporary scalars
76      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! temporary scalars
77      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu, zfv
78      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw
79      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zlu_uu, zlv_vv, zlu_uv, zlv_vu
80      !!----------------------------------------------------------------------
81      !
82      IF( nn_timing == 1 )  CALL timing_start('dyn_adv_ubs')
83      !
84      CALL wrk_alloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
85      CALL wrk_alloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                               )
86      !
87      IF( kt == nit000 ) THEN
88         IF(lwp) WRITE(numout,*)
89         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
90         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
91      ENDIF
92      !
93      zfu_t(:,:,:) = 0._wp
94      zfv_t(:,:,:) = 0._wp
95      zfu_f(:,:,:) = 0._wp
96      zfv_f(:,:,:) = 0._wp
97      !
98      zlu_uu(:,:,:,:) = 0._wp
99      zlv_vv(:,:,:,:) = 0._wp 
100      zlu_uv(:,:,:,:) = 0._wp 
101      zlv_vu(:,:,:,:) = 0._wp 
102      !
103      IF( l_trddyn ) THEN           ! trends: store the input trends
104         zfu_uw(:,:,:) = ua(:,:,:)
105         zfv_vw(:,:,:) = va(:,:,:)
106      ENDIF
107      !                                      ! =========================== !
108      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
109         !                                   ! =========================== !
110         !                                         ! horizontal volume fluxes
111         zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
112         zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
113         !           
114         DO jj = 2, jpjm1                          ! laplacian
115            DO ji = fs_2, fs_jpim1   ! vector opt.
116               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)
117               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)
118               zlu_uv(ji,jj,jk,1) = ( ub (ji  ,jj+1,jk) - ub (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
119                  &               - ( ub (ji  ,jj  ,jk) - ub (ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
120               zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj  ,jk) - vb (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
121                  &               - ( vb (ji  ,jj  ,jk) - vb (ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
122               !
123               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)
124               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)
125               zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
126                  &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
127               zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
128                  &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
129            END DO
130         END DO
131      END DO
132      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
133      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
134      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
135      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 
136      !
137      !                                      ! ====================== !
138      !                                      !  Horizontal advection  !
139      DO jk = 1, jpkm1                       ! ====================== !
140         !                                         ! horizontal volume fluxes
141         zfu(:,:,jk) = 0.25 * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
142         zfv(:,:,jk) = 0.25 * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
143         !
144         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
145            DO ji = 1, fs_jpim1   ! vector opt.
146               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
147               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
148               !
149               IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
150               ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
151               ENDIF
152               IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
153               ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
154               ENDIF
155               !
156               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
157                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
158                  &                * ( zui - gamma1 * zl_u)
159               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
160                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
161                  &                * ( zvj - gamma1 * zl_v)
162               !
163               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
164               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
165               IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
166               ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
167               ENDIF
168               IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
169               ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
170               ENDIF
171               !
172               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
173                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
174               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
175                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
176            END DO
177         END DO
178         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
179            DO ji = fs_2, fs_jpim1   ! vector opt.
180               zbu = e1e2u(ji,jj) * e3u_n(ji,jj,jk)
181               zbv = e1e2v(ji,jj) * e3v_n(ji,jj,jk)
182               !
183               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
184                  &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
185               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    &
186                  &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
187            END DO
188         END DO
189      END DO
190      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
191         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
192         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
193         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
194         zfu_t(:,:,:) = ua(:,:,:)
195         zfv_t(:,:,:) = va(:,:,:)
196      ENDIF
197      !                                      ! ==================== !
198      !                                      !  Vertical advection  !
199      !                                      ! ==================== !
200      DO jj = 2, jpjm1                             ! surface/bottom advective fluxes set to zero                 
201         DO ji = fs_2, fs_jpim1
202            zfu_uw(ji,jj,jpk) = 0._wp
203            zfv_vw(ji,jj,jpk) = 0._wp
204            zfu_uw(ji,jj, 1 ) = 0._wp
205            zfv_vw(ji,jj, 1 ) = 0._wp
206         END DO
207      END DO
208      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
209         DO jj = 2, jpjm1
210            DO ji = fs_2, fs_jpim1
211               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)
212               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)
213            END DO
214         END DO
215      ENDIF
216      DO jk = 2, jpkm1                          ! interior fluxes
217         DO jj = 2, jpjm1
218            DO ji = fs_2, fs_jpim1
219               zfw(ji,jj,jk) = 0.25 * e1e2t(ji,jj) * wn(ji,jj,jk)
220            END DO
221         END DO
222         DO jj = 2, jpjm1
223            DO ji = fs_2, fs_jpim1   ! vector opt.
224               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
225               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
226            END DO
227         END DO
228      END DO
229      DO jk = 1, jpkm1                          ! divergence of vertical momentum flux divergence
230         DO jj = 2, jpjm1
231            DO ji = fs_2, fs_jpim1   ! vector opt.
232               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)
233               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)
234            END DO
235         END DO
236      END DO
237      !
238      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
239         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
240         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
241         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
242      ENDIF
243      !                                         ! Control print
244      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
245         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
246      !
247      CALL wrk_dealloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
248      CALL wrk_dealloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                               )
249      !
250      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_ubs')
251      !
252   END SUBROUTINE dyn_adv_ubs
253
254   !!==============================================================================
255END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.