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

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by gm, 12 years ago

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