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

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

#1613: vvl by default: add ln_linssh and remove key_vvl
  • Property svn:keywords set to Id
File size: 13.9 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           ! Save ua and va trends
104         zfu_uw(:,:,:) = ua(:,:,:)
105         zfv_vw(:,:,:) = va(:,:,:)
106      ENDIF
107
108      !                                      ! =========================== !
109      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
110         !                                   ! =========================== !
111         !                                         ! horizontal volume fluxes
112         zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
113         zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
114         !           
115         DO jj = 2, jpjm1                          ! laplacian
116            DO ji = fs_2, fs_jpim1   ! vector opt.
117               !
118               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)
119               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)
120               zlu_uv(ji,jj,jk,1) = ( ub (ji  ,jj+1,jk) - ub (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
121                  &               - ( ub (ji  ,jj  ,jk) - ub (ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
122               zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj  ,jk) - vb (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
123                  &               - ( vb (ji  ,jj  ,jk) - vb (ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
124               !
125               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)
126               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)
127               zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
128                  &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
129               zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
130                  &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
131            END DO
132         END DO
133      END DO
134      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
135      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
136      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
137      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 
138     
139      !                                      ! ====================== !
140      !                                      !  Horizontal advection  !
141      DO jk = 1, jpkm1                       ! ====================== !
142         !                                         ! horizontal volume fluxes
143         zfu(:,:,jk) = 0.25 * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
144         zfv(:,:,jk) = 0.25 * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
145         !
146         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
147            DO ji = 1, fs_jpim1   ! vector opt.
148               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
149               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
150               !
151               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
152               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
153               ENDIF
154               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
155               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
156               ENDIF
157               !
158               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
159                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
160                  &                * ( zui - gamma1 * zl_u)
161               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
162                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
163                  &                * ( zvj - gamma1 * zl_v)
164               !
165               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
166               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
167               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
168               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
169               ENDIF
170               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
171               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
172               ENDIF
173               !
174               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
175                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
176               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
177                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
178            END DO
179         END DO
180         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
181            DO ji = fs_2, fs_jpim1   ! vector opt.
182               zbu = e1e2u(ji,jj) * e3u_n(ji,jj,jk)
183               zbv = e1e2v(ji,jj) * e3v_n(ji,jj,jk)
184               !
185               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
186                  &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
187               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    &
188                  &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
189            END DO
190         END DO
191      END DO
192      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
193         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
194         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
195         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
196         zfu_t(:,:,:) = ua(:,:,:)
197         zfv_t(:,:,:) = va(:,:,:)
198      ENDIF
199
200      !                                      ! ==================== !
201      !                                      !  Vertical advection  !
202      DO jk = 1, jpkm1                       ! ==================== !
203         !                                         ! Vertical volume fluxesÊ
204         zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk)
205         !
206         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
207            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
208            zfv_vw(:,:,jpk) = 0.e0
209            !                                           ! Surface value :
210            IF( .NOT.ln_linssh ) THEN                        ! variable volume : flux set to zero
211               zfu_uw(:,:, 1 ) = 0._wp
212               zfv_vw(:,:, 1 ) = 0._wp
213            ELSE                                             ! constant volume : advection through the surface
214               DO jj = 2, jpjm1
215                  DO ji = fs_2, fs_jpim1
216                     zfu_uw(ji,jj, 1 ) = 2._wp * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
217                     zfv_vw(ji,jj, 1 ) = 2._wp * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
218                  END DO
219               END DO
220            ENDIF
221         ELSE                                      ! interior fluxes
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         ENDIF
229      END DO
230      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
231         DO jj = 2, jpjm1 
232            DO ji = fs_2, fs_jpim1   ! vector opt.
233               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
234                  &  / ( e1e2u(ji,jj) * e3u_n(ji,jj,jk) )
235               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
236                  &  / ( e1e2v(ji,jj) * e3v_n(ji,jj,jk) )
237            END DO
238         END DO
239      END DO
240      !
241      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
242         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
243         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
244         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
245      ENDIF
246      !                                            ! Control print
247      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
248         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
249      !
250      CALL wrk_dealloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
251      CALL wrk_dealloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                               )
252      !
253      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_ubs')
254      !
255   END SUBROUTINE dyn_adv_ubs
256
257   !!==============================================================================
258END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.