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/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 5038

Last change on this file since 5038 was 5038, checked in by jamesharle, 9 years ago

Merging branch with HEAD of the trunk

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