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

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

Last change on this file since 5062 was 5062, checked in by flavoni, 9 years ago

correction to previous committ
  • Property svn:keywords set to Id
File size: 14.2 KB
RevLine 
[643]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   !!======================================================================
[1566]7   !! History :  2.0  ! 2006-08  (R. Benshila, L. Debreu)  Original code
8   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
[643]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
[5062]18   USE trd_oce        ! trends: ocean variables
19   USE trddyn         ! trend manager: dynamics
20   !
[2715]21   USE in_out_manager ! I/O manager
[1129]22   USE prtctl         ! Print control
[2715]23   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
24   USE lib_mpp        ! MPP library
[5062]25   USE wrk_nemo       ! Memory Allocation
26   USE timing         ! Timing
[643]27
28   IMPLICIT NONE
29   PRIVATE
30
[3294]31   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS
[4153]32   REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0   2nd order  ; =1/32 4th order centred
[643]33
[1566]34   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
[643]35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
[2715]40   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[1152]41   !! $Id$
[2715]42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[643]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
[1566]51      !!              and the general trend of the momentum equation.
[643]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
[4153]60      !!                       = 1/32 4th order finite differencing
[643]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
[4153]67      !!      gamma1=1/3 and gamma2=1/32.
[643]68      !!
[1566]69      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
[643]70      !!
71      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
72      !!----------------------------------------------------------------------
73      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
[2715]74      !
[1566]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
[3294]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
[643]81      !!----------------------------------------------------------------------
[3294]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      !
[643]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
[3294]93      !
[2715]94      zfu_t(:,:,:) = 0._wp
95      zfv_t(:,:,:) = 0._wp
96      zfu_f(:,:,:) = 0._wp
97      zfv_f(:,:,:) = 0._wp
[1566]98      !
[2715]99      zlu_uu(:,:,:,:) = 0._wp
100      zlv_vv(:,:,:,:) = 0._wp 
101      zlu_uv(:,:,:,:) = 0._wp 
102      zlv_vu(:,:,:,:) = 0._wp 
[643]103
[1129]104      IF( l_trddyn ) THEN           ! Save ua and va trends
105         zfu_uw(:,:,:) = ua(:,:,:)
106         zfv_vw(:,:,:) = va(:,:,:)
107      ENDIF
108
[1566]109      !                                      ! =========================== !
110      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
111         !                                   ! =========================== !
112         !                                         ! horizontal volume fluxes
[643]113         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
114         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]115         !           
116         DO jj = 2, jpjm1                          ! laplacian
[643]117            DO ji = fs_2, fs_jpim1   ! vector opt.
[5061]118               !
[643]119               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)
[5061]120               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)
121               zlu_uv(ji,jj,jk,1) = ( ub (ji,jj+1,jk)-ub (ji,jj,jk) ) * fmask(ji,jj,jk)-( 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)-( vb (ji,jj,jk)-vb (ji-1,jj,jk) ) * fmask(ji-1,jj,jk)
[2715]123               !
[643]124               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)
125               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)
[5061]126               zlu_uv(ji,jj,jk,2) = ( zfu (ji,jj+1,jk)-zfu (ji,jj,jk) ) * fmask(ji,jj,jk)-( 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)-( zfv (ji,jj,jk)-zfv (ji-1,jj,jk) ) * fmask(ji-1,jj,jk)
[643]128            END DO
129         END DO
[1566]130      END DO
131!!gm BUG !!!  just below this should be +1 in all the communications
[2715]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!
[1566]137!!gm corrected:
138      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
139      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
140      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
141      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 
142!!gm end
143     
144      !                                      ! ====================== !
145      !                                      !  Horizontal advection  !
146      DO jk = 1, jpkm1                       ! ====================== !
147         !                                         ! horizontal volume fluxes
[643]148         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
149         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]150         !
151         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
[643]152            DO ji = 1, fs_jpim1   ! vector opt.
153               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
154               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
[1566]155               !
[643]156               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
157               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
158               ENDIF
159               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
160               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
161               ENDIF
[1566]162               !
[643]163               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
164                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
165                  &                * ( zui - gamma1 * zl_u)
166               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
167                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
168                  &                * ( zvj - gamma1 * zl_v)
[1566]169               !
[643]170               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
171               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
172               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
173               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
174               ENDIF
175               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
176               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
177               ENDIF
[1566]178               !
[643]179               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
180                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
181               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
182                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
183            END DO
184         END DO
[1566]185         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
[643]186            DO ji = fs_2, fs_jpim1   ! vector opt.
187               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
188               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
[1566]189               !
190               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
191                  &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
192               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    &
193                  &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
[643]194            END DO
195         END DO
[1566]196      END DO
197      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
[1129]198         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
199         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
[5062]200         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
[1129]201         zfu_t(:,:,:) = ua(:,:,:)
202         zfv_t(:,:,:) = va(:,:,:)
203      ENDIF
204
[1566]205      !                                      ! ==================== !
206      !                                      !  Vertical advection  !
207      DO jk = 1, jpkm1                       ! ==================== !
208         !                                         ! Vertical volume fluxesÊ
[643]209         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
[1566]210         !
211         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
212            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
[643]213            zfv_vw(:,:,jpk) = 0.e0
[1566]214            !                                           ! Surface value :
215            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
[643]216               zfu_uw(:,:, 1 ) = 0.e0   
217               zfv_vw(:,:, 1 ) = 0.e0
[1566]218            ELSE                                             ! constant volume : advection through the surface
[643]219               DO jj = 2, jpjm1
220                  DO ji = fs_2, fs_jpim1
221                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
222                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
223                  END DO
224               END DO
225            ENDIF
[1566]226         ELSE                                      ! interior fluxes
[643]227            DO jj = 2, jpjm1
228               DO ji = fs_2, fs_jpim1   ! vector opt.
229                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
230                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
231               END DO
232            END DO
233         ENDIF
234      END DO
[1566]235      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
[643]236         DO jj = 2, jpjm1 
237            DO ji = fs_2, fs_jpim1   ! vector opt.
[1566]238               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
[643]239                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
[1566]240               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
[643]241                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
242            END DO
243         END DO
244      END DO
[1566]245      !
246      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
[1129]247         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
248         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
[5062]249         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
[1129]250      ENDIF
[1566]251      !                                            ! Control print
[1129]252      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
253         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]254      !
[3294]255      CALL wrk_dealloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
256      CALL wrk_dealloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                               )
[2715]257      !
[3294]258      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_ubs')
259      !
[643]260   END SUBROUTINE dyn_adv_ubs
261
262   !!==============================================================================
263END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.