source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 6750

Last change on this file since 6750 was 6750, checked in by jchanut, 4 years ago

Reproducibility issue with flux form advection, ticket #1749

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