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 @ 4460

Last change on this file since 4460 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
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
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
31   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
32
33   !! * Control permutation of array indices
34#  include "oce_ftrans.h90"
35#  include "dom_oce_ftrans.h90"
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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
52      !!              and the general trend of the momentum equation.
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      !!
70      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
71      !!
72      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
73      !!----------------------------------------------------------------------
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
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
95      !
96      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
97      !
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
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
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
118      !
119      zlu_uu(:,:,:,:) = 0._wp
120      zlv_vv(:,:,:,:) = 0._wp 
121      zlu_uv(:,:,:,:) = 0._wp 
122      zlv_vu(:,:,:,:) = 0._wp 
123
124      IF( l_trddyn ) THEN           ! Save ua and va trends
125         zfu_uw(:,:,:) = ua(:,:,:)
126         zfv_vw(:,:,:) = va(:,:,:)
127      ENDIF
128
129      !                                      ! =========================== !
130      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
131         !                                   ! =========================== !
132         !                                         ! horizontal volume fluxes
133         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
134         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
135         !           
136         DO jj = 2, jpjm1                          ! laplacian
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)
142               !
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
149      END DO
150!!gm BUG !!!  just below this should be +1 in all the communications
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!
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
167         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
168         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
169         !
170         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
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) )
174               !
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
181               !
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)
188               !
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
197               !
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
204         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
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)
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
213            END DO
214         END DO
215      END DO
216      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
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
224      !                                      ! ==================== !
225      !                                      !  Vertical advection  !
226      DO jk = 1, jpkm1                       ! ==================== !
227         !                                         ! Vertical volume fluxesÊ
228         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
229         !
230         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
231            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
232            zfv_vw(:,:,jpk) = 0.e0
233            !                                           ! Surface value :
234            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
235               zfu_uw(:,:, 1 ) = 0.e0   
236               zfv_vw(:,:, 1 ) = 0.e0
237            ELSE                                             ! constant volume : advection through the surface
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
245         ELSE                                      ! interior fluxes
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
254      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
255         DO jj = 2, jpjm1 
256            DO ji = fs_2, fs_jpim1   ! vector opt.
257               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
258                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
259               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
260                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
261            END DO
262         END DO
263      END DO
264      !
265      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
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
270      !                                            ! Control print
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' )
273      !
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      !
277   END SUBROUTINE dyn_adv_ubs
278
279   !!==============================================================================
280END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.