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/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 3186

Last change on this file since 3186 was 3186, checked in by smasson, 12 years ago

dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo
  • Property svn:keywords set to Id
File size: 14.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 trdmod         ! ocean dynamics trends
19   USE trdmod_oce     ! ocean variables trends
20   USE in_out_manager ! I/O manager
21   USE prtctl         ! Print control
22   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
23   USE lib_mpp        ! MPP library
24   USE wrk_nemo        ! Memory Allocation
25   USE timing          ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   REAL(wp), PARAMETER :: gamma1 = 1._wp/4._wp  ! =1/4 quick      ; =1/3  3rd order UBS
31   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred
32
33   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
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/8  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/4 and gamma2=1/8.
67      !!
68      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
69      !!
70      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
71      !!----------------------------------------------------------------------
72      USE oce     , ONLY:   tsa             ! tsa used as 2 3D workspace
73      !
74      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
75      !
76      INTEGER  ::   ji, jj, jk            ! dummy loop indices
77      REAL(wp) ::   zbu, zbv    ! temporary scalars
78      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! temporary scalars
79      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu, zfv
80      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw
81      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zlu_uu, zlv_vv, zlu_uv, zlv_vu
82      !!----------------------------------------------------------------------
83      !
84      IF( nn_timing == 1 )  CALL timing_start('dyn_adv_ubs')
85      !
86      CALL wrk_alloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfw )
87      CALL wrk_alloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                     )
88      zfu => tsa(:,:,:,1) 
89      zfv => tsa(:,:,:,2) 
90      !
91      IF( kt == nit000 ) THEN
92         IF(lwp) WRITE(numout,*)
93         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
94         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
95      ENDIF
96      !
97      zfu_t(:,:,:) = 0._wp
98      zfv_t(:,:,:) = 0._wp
99      zfu_f(:,:,:) = 0._wp
100      zfv_f(:,:,:) = 0._wp
101      !
102      zlu_uu(:,:,:,:) = 0._wp
103      zlv_vv(:,:,:,:) = 0._wp 
104      zlu_uv(:,:,:,:) = 0._wp 
105      zlv_vu(:,:,:,:) = 0._wp 
106
107      IF( l_trddyn ) THEN           ! Save ua and va trends
108         zfu_uw(:,:,:) = ua(:,:,:)
109         zfv_vw(:,:,:) = va(:,:,:)
110      ENDIF
111
112      !                                      ! =========================== !
113      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
114         !                                   ! =========================== !
115         !                                         ! horizontal volume fluxes
116         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
117         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
118         !           
119         DO jj = 2, jpjm1                          ! laplacian
120            DO ji = fs_2, fs_jpim1   ! vector opt.
121               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)
122               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) 
123               zlu_uv(ji,jj,jk,1) = ( ub (ji,jj+1,jk)-2.*ub (ji,jj,jk)+ub (ji,jj-1,jk) ) * umask(ji,jj,jk)
124               zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj,jk)-2.*vb (ji,jj,jk)+vb (ji-1,jj,jk) ) * vmask(ji,jj,jk)
125               !
126               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)
127               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)
128               zlu_uv(ji,jj,jk,2) = ( zfu(ji,jj+1,jk)-2.*zfu(ji,jj,jk)+zfu(ji,jj-1,jk) ) * umask(ji,jj,jk)
129               zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj,jk)-2.*zfv(ji,jj,jk)+zfv(ji-1,jj,jk) ) * vmask(ji,jj,jk)
130            END DO
131         END DO
132      END DO
133!!gm BUG !!!  just below this should be +1 in all the communications
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!!gm corrected:
140      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
141      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
142      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
143      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 
144!!gm end
145     
146      !                                      ! ====================== !
147      !                                      !  Horizontal advection  !
148      DO jk = 1, jpkm1                       ! ====================== !
149         !                                         ! horizontal volume fluxes
150         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
151         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
152         !
153         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
154            DO ji = 1, fs_jpim1   ! vector opt.
155               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
156               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
157               !
158               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
159               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
160               ENDIF
161               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
162               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
163               ENDIF
164               !
165               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
166                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
167                  &                * ( zui - gamma1 * zl_u)
168               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
169                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
170                  &                * ( zvj - gamma1 * zl_v)
171               !
172               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
173               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
174               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
175               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
176               ENDIF
177               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
178               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
179               ENDIF
180               !
181               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
182                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
183               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
184                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
185            END DO
186         END DO
187         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
188            DO ji = fs_2, fs_jpim1   ! vector opt.
189               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
190               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
191               !
192               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
193                  &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
194               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    &
195                  &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
196            END DO
197         END DO
198      END DO
199      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
200         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
201         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
202         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
203         zfu_t(:,:,:) = ua(:,:,:)
204         zfv_t(:,:,:) = va(:,:,:)
205      ENDIF
206
207      !                                      ! ==================== !
208      !                                      !  Vertical advection  !
209      DO jk = 1, jpkm1                       ! ==================== !
210         !                                         ! Vertical volume fluxesÊ
211         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
212         !
213         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
214            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
215            zfv_vw(:,:,jpk) = 0.e0
216            !                                           ! Surface value :
217            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
218               zfu_uw(:,:, 1 ) = 0.e0   
219               zfv_vw(:,:, 1 ) = 0.e0
220            ELSE                                             ! constant volume : advection through the surface
221               DO jj = 2, jpjm1
222                  DO ji = fs_2, fs_jpim1
223                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
224                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
225                  END DO
226               END DO
227            ENDIF
228         ELSE                                      ! interior fluxes
229            DO jj = 2, jpjm1
230               DO ji = fs_2, fs_jpim1   ! vector opt.
231                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
232                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
233               END DO
234            END DO
235         ENDIF
236      END DO
237      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
238         DO jj = 2, jpjm1 
239            DO ji = fs_2, fs_jpim1   ! vector opt.
240               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
241                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
242               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
243                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
244            END DO
245         END DO
246      END DO
247      !
248      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
249         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
250         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
251         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
252      ENDIF
253      !                                            ! Control print
254      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
255         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
256      !
257      CALL wrk_dealloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfw )
258      CALL wrk_dealloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                     )
259      !
260      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_ubs')
261      !
262   END SUBROUTINE dyn_adv_ubs
263
264   !!==============================================================================
265END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.