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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 4400

Last change on this file since 4400 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 14.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
[1129]18   USE trdmod         ! ocean dynamics trends
19   USE trdmod_oce     ! ocean variables trends
[2715]20   USE in_out_manager ! I/O manager
[1129]21   USE prtctl         ! Print control
[2715]22   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
23   USE lib_mpp        ! MPP library
[643]24
25   IMPLICIT NONE
26   PRIVATE
27
28   REAL(wp), PARAMETER :: gamma1 = 1._wp/4._wp  ! =1/4 quick      ; =1/3  3rd order UBS
29   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred
30
[1566]31   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
[643]32
[3211]33   !! * Control permutation of array indices
34#  include "oce_ftrans.h90"
35#  include "dom_oce_ftrans.h90"
36
[643]37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
[2715]41   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[1152]42   !! $Id$
[2715]43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[643]44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE dyn_adv_ubs( kt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE dyn_adv_ubs  ***
50      !!
51      !! ** Purpose :   Compute the now momentum advection trend in flux form
[1566]52      !!              and the general trend of the momentum equation.
[643]53      !!
54      !! ** Method  :   The scheme is the one implemeted in ROMS. It depends
55      !!      on two parameter gamma1 and gamma2. The former control the
56      !!      upstream baised part of the scheme and the later the centred
57      !!      part:     gamma1 = 0    pure centered  (no diffusive part)
58      !!                       = 1/4  Quick scheme
59      !!                       = 1/3  3rd order Upstream biased scheme
60      !!                gamma2 = 0    2nd order finite differencing
61      !!                       = 1/8  4th order finite differencing
62      !!      For stability reasons, the first term of the fluxes which cor-
63      !!      responds to a second order centered scheme is evaluated using 
64      !!      the now velocity (centered in time) while the second term which 
65      !!      is the diffusive part of the scheme, is evaluated using the
66      !!      before velocity (forward in time).
67      !!      Default value (hard coded in the begining of the module) are
68      !!      gamma1=1/4 and gamma2=1/8.
69      !!
[1566]70      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
[643]71      !!
72      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
73      !!----------------------------------------------------------------------
[2715]74      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
75      USE oce     , ONLY:   zfu    => ta       , zfv    => sa      ! (ta,sa) used as 3D workspace
76      USE wrk_nemo, ONLY:   zfu_t  => wrk_3d_1 , zfv_t  =>wrk_3d_4 , zfu_uw =>wrk_3d_6   ! 3D workspace
77      USE wrk_nemo, ONLY:   zfu_f  => wrk_3d_2 , zfv_f  =>wrk_3d_5 , zfv_vw =>wrk_3d_7
78      USE wrk_nemo, ONLY:   zfw    => wrk_3d_3
79      USE wrk_nemo, ONLY:   zlu_uu => wrk_4d_1 , zlv_vv=>wrk_4d_3   ! 4D workspace
80      USE wrk_nemo, ONLY:   zlu_uv => wrk_4d_2 , zlv_vu=>wrk_4d_4
[3211]81      !! DCSE_NEMO: module variables renamed, need additional directives
82!FTRANS zfu :I :I :z
83!FTRANS zfv :I :I :z
84!FTRANS zfu_t :I :I :z
85!FTRANS zfv_t :I :I :z
86!FTRANS zfu_uw :I :I :z
87!FTRANS zfu_f :I :I :z
88!FTRANS zfv_f :I :I :z
89!FTRANS zfv_vw :I :I :z
90!FTRANS zfw :I :I :z
91!FTRANS zlu_uu :I :I :z :I
92!FTRANS zlv_vv :I :I :z :I
93!FTRANS zlu_uv :I :I :z :I
94!FTRANS zlv_vu :I :I :z :I
[2715]95      !
[643]96      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
[2715]97      !
[1566]98      INTEGER  ::   ji, jj, jk            ! dummy loop indices
99      REAL(wp) ::   zbu, zbv    ! temporary scalars
100      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! temporary scalars
[643]101      !!----------------------------------------------------------------------
102
103      IF( kt == nit000 ) THEN
104         IF(lwp) WRITE(numout,*)
105         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
106         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
107      ENDIF
[2715]108
109      ! Check that required workspace arrays are not already in use
110      IF( wrk_in_use(3, 1,2,3,4,5,6,7) .OR. wrk_in_use(4, 1,2,3,4) ) THEN
111         CALL ctl_stop('dyn_adv_ubs: requested workspace array unavailable')   ;   RETURN
112      ENDIF
113
114      zfu_t(:,:,:) = 0._wp
115      zfv_t(:,:,:) = 0._wp
116      zfu_f(:,:,:) = 0._wp
117      zfv_f(:,:,:) = 0._wp
[1566]118      !
[2715]119      zlu_uu(:,:,:,:) = 0._wp
120      zlv_vv(:,:,:,:) = 0._wp 
121      zlu_uv(:,:,:,:) = 0._wp 
122      zlv_vu(:,:,:,:) = 0._wp 
[643]123
[1129]124      IF( l_trddyn ) THEN           ! Save ua and va trends
125         zfu_uw(:,:,:) = ua(:,:,:)
126         zfv_vw(:,:,:) = va(:,:,:)
127      ENDIF
128
[1566]129      !                                      ! =========================== !
130      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
131         !                                   ! =========================== !
132         !                                         ! horizontal volume fluxes
[643]133         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
134         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]135         !           
136         DO jj = 2, jpjm1                          ! laplacian
[643]137            DO ji = fs_2, fs_jpim1   ! vector opt.
138               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)
139               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) 
140               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)
141               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)
[2715]142               !
[643]143               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)
144               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)
145               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)
146               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)
147            END DO
148         END DO
[1566]149      END DO
150!!gm BUG !!!  just below this should be +1 in all the communications
[2715]151!      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.)
152!      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.)
153!      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.)
154!      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.)
155!
[1566]156!!gm corrected:
157      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
158      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
159      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
160      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 
161!!gm end
162     
163      !                                      ! ====================== !
164      !                                      !  Horizontal advection  !
165      DO jk = 1, jpkm1                       ! ====================== !
166         !                                         ! horizontal volume fluxes
[643]167         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
168         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]169         !
170         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
[643]171            DO ji = 1, fs_jpim1   ! vector opt.
172               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
173               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
[1566]174               !
[643]175               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
176               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
177               ENDIF
178               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
179               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
180               ENDIF
[1566]181               !
[643]182               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
183                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
184                  &                * ( zui - gamma1 * zl_u)
185               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
186                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
187                  &                * ( zvj - gamma1 * zl_v)
[1566]188               !
[643]189               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
190               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
191               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
192               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
193               ENDIF
194               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
195               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
196               ENDIF
[1566]197               !
[643]198               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
199                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
200               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
201                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
202            END DO
203         END DO
[1566]204         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
[643]205            DO ji = fs_2, fs_jpim1   ! vector opt.
206               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
207               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
[1566]208               !
209               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
210                  &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
211               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    &
212                  &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
[643]213            END DO
214         END DO
[1566]215      END DO
216      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
[1129]217         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
218         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
219         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
220         zfu_t(:,:,:) = ua(:,:,:)
221         zfv_t(:,:,:) = va(:,:,:)
222      ENDIF
223
[1566]224      !                                      ! ==================== !
225      !                                      !  Vertical advection  !
226      DO jk = 1, jpkm1                       ! ==================== !
227         !                                         ! Vertical volume fluxesÊ
[643]228         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
[1566]229         !
230         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
231            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
[643]232            zfv_vw(:,:,jpk) = 0.e0
[1566]233            !                                           ! Surface value :
234            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
[643]235               zfu_uw(:,:, 1 ) = 0.e0   
236               zfv_vw(:,:, 1 ) = 0.e0
[1566]237            ELSE                                             ! constant volume : advection through the surface
[643]238               DO jj = 2, jpjm1
239                  DO ji = fs_2, fs_jpim1
240                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
241                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
242                  END DO
243               END DO
244            ENDIF
[1566]245         ELSE                                      ! interior fluxes
[643]246            DO jj = 2, jpjm1
247               DO ji = fs_2, fs_jpim1   ! vector opt.
248                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
249                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
250               END DO
251            END DO
252         ENDIF
253      END DO
[1566]254      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
[643]255         DO jj = 2, jpjm1 
256            DO ji = fs_2, fs_jpim1   ! vector opt.
[1566]257               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
[643]258                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
[1566]259               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
[643]260                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
261            END DO
262         END DO
263      END DO
[1566]264      !
265      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
[1129]266         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
267         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
268         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
269      ENDIF
[1566]270      !                                            ! Control print
[1129]271      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
272         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]273      !
[2715]274      IF( wrk_not_released(3, 1,2,3,4,5,6,7) .OR.   &
275          wrk_not_released(4, 1,2,3,4)       )   CALL ctl_stop('dyn_adv_ubs: failed to release workspace array')
276      !
[643]277   END SUBROUTINE dyn_adv_ubs
278
279   !!==============================================================================
280END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.