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

source: trunk/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 1950

Last change on this file since 1950 was 1566, checked in by rblod, 15 years ago

Cosmetic changes: suppress useless variables and code review of the code changed when suppressing rigid-lid, see ticket #508
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Id
File size: 13.6 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
18   USE in_out_manager ! I/O manager
19   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
[1129]20   USE trdmod         ! ocean dynamics trends
21   USE trdmod_oce     ! ocean variables trends
22   USE prtctl         ! Print control
[643]23
24   IMPLICIT NONE
25   PRIVATE
26
27   REAL(wp), PARAMETER :: gamma1 = 1._wp/4._wp  ! =1/4 quick      ; =1/3  3rd order UBS
28   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred
29
[1566]30   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
[643]31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
[1566]36   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
[1152]37   !! $Id$
[643]38   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE dyn_adv_ubs( kt )
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE dyn_adv_ubs  ***
46      !!
47      !! ** Purpose :   Compute the now momentum advection trend in flux form
[1566]48      !!              and the general trend of the momentum equation.
[643]49      !!
50      !! ** Method  :   The scheme is the one implemeted in ROMS. It depends
51      !!      on two parameter gamma1 and gamma2. The former control the
52      !!      upstream baised part of the scheme and the later the centred
53      !!      part:     gamma1 = 0    pure centered  (no diffusive part)
54      !!                       = 1/4  Quick scheme
55      !!                       = 1/3  3rd order Upstream biased scheme
56      !!                gamma2 = 0    2nd order finite differencing
57      !!                       = 1/8  4th order finite differencing
58      !!      For stability reasons, the first term of the fluxes which cor-
59      !!      responds to a second order centered scheme is evaluated using 
60      !!      the now velocity (centered in time) while the second term which 
61      !!      is the diffusive part of the scheme, is evaluated using the
62      !!      before velocity (forward in time).
63      !!      Default value (hard coded in the begining of the module) are
64      !!      gamma1=1/4 and gamma2=1/8.
65      !!
[1566]66      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
[643]67      !!
68      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
69      !!----------------------------------------------------------------------
[1566]70      USE oce, ONLY:   zfu => ta   ! use ta as 3D workspace
71      USE oce, ONLY:   zfv => sa   ! use sa as 3D workspace
72      !!
[643]73      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
[1566]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
[643]78      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f     ! temporary workspace
79      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f     !    "           "
80      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfw, zfu_uw, zfv_vw
81      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv   ! temporary workspace
82      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu   ! temporary workspace
83      !!----------------------------------------------------------------------
84
85      IF( kt == nit000 ) THEN
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
88         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
89      ENDIF
90      zfu_t(:,:,:) = 0.e0
91      zfv_t(:,:,:) = 0.e0
92      zfu_f(:,:,:) = 0.e0
93      zfv_f(:,:,:) = 0.e0
[1566]94      !
[643]95      zlu_uu(:,:,:,:) = 0.e0 
96      zlv_vv(:,:,:,:) = 0.e0 
97      zlu_uv(:,:,:,:) = 0.e0 
98      zlv_vu(:,:,:,:) = 0.e0 
99
[1129]100      IF( l_trddyn ) THEN           ! Save ua and va trends
101         zfu_uw(:,:,:) = ua(:,:,:)
102         zfv_vw(:,:,:) = va(:,:,:)
103      ENDIF
104
[1566]105      !                                      ! =========================== !
106      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
107         !                                   ! =========================== !
108         !                                         ! horizontal volume fluxes
[643]109         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
110         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]111         !           
112         DO jj = 2, jpjm1                          ! laplacian
[643]113            DO ji = fs_2, fs_jpim1   ! vector opt.
114               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)
115               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) 
116               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)
117               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)
118               
119               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)
120               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)
121               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)
122               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)
123            END DO
124         END DO
[1566]125      END DO
126!!gm BUG !!!  just below this should be +1 in all the communications
127      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.)
128      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.)
129      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.)
130      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 
[643]131
[1566]132!!gm corrected:
133      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
134      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
135      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
136      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 
137!!gm end
138     
139      !                                      ! ====================== !
140      !                                      !  Horizontal advection  !
141      DO jk = 1, jpkm1                       ! ====================== !
142         !                                         ! horizontal volume fluxes
[643]143         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
144         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]145         !
146         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
[643]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) )
[1566]150               !
[643]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
[1566]157               !
[643]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)
[1566]164               !
[643]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
[1566]173               !
[643]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
[1566]180         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
[643]181            DO ji = fs_2, fs_jpim1   ! vector opt.
182               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
183               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
[1566]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
[643]189            END DO
190         END DO
[1566]191      END DO
192      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
[1129]193         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
194         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
195         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
196         zfu_t(:,:,:) = ua(:,:,:)
197         zfv_t(:,:,:) = va(:,:,:)
198      ENDIF
199
[1566]200      !                                      ! ==================== !
201      !                                      !  Vertical advection  !
202      DO jk = 1, jpkm1                       ! ==================== !
203         !                                         ! Vertical volume fluxesÊ
[643]204         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
[1566]205         !
206         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
207            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
[643]208            zfv_vw(:,:,jpk) = 0.e0
[1566]209            !                                           ! Surface value :
210            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
[643]211               zfu_uw(:,:, 1 ) = 0.e0   
212               zfv_vw(:,:, 1 ) = 0.e0
[1566]213            ELSE                                             ! constant volume : advection through the surface
[643]214               DO jj = 2, jpjm1
215                  DO ji = fs_2, fs_jpim1
216                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
217                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
218                  END DO
219               END DO
220            ENDIF
[1566]221         ELSE                                      ! interior fluxes
[643]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
[1566]230      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
[643]231         DO jj = 2, jpjm1 
232            DO ji = fs_2, fs_jpim1   ! vector opt.
[1566]233               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
[643]234                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
[1566]235               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
[643]236                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
237            END DO
238         END DO
239      END DO
[1566]240      !
241      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
[1129]242         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
243         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
244         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
245      ENDIF
[1566]246      !                                            ! Control print
[1129]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' )
[1566]249      !
[643]250   END SUBROUTINE dyn_adv_ubs
251
252   !!==============================================================================
253END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.