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

Last change on this file since 1129 was 1129, checked in by ctlod, 16 years ago

trunk: compilation error when key_trddyn is active related to the dynamics advection flux formulation, see ticket: #211
  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 13.5 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 :  9.0  !  06-08  (R. Benshila, L. Debreu)  Original code
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   dyn_adv_ubs   : flux form momentum advection using    (ln_dynadv=T)
12   !!                   an 3rd order Upstream Biased Scheme or Quick scheme
13   !!                   combined with 2nd or 4th order finite differences
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers
16   USE dom_oce        ! ocean space and time domain
17   USE dynspg_oce     ! surface pressure gradient
18   USE in_out_manager ! I/O manager
19   USE dynspg_rl      ! surface pressure gradient
20   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
21   USE trdmod         ! ocean dynamics trends
22   USE trdmod_oce     ! ocean variables trends
23   USE prtctl         ! Print control
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   !! * Routine accessibility
32   PUBLIC dyn_adv_ubs                            ! routine called by step.F90
33
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !!   OPA 9.0 , LODYC-IPSL  (2006)
39   !! $Header$
40   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
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/4 and gamma2=1/8.
67      !!
68      !! ** Action : - update (ua,va) with the 3D advective momentum trends
69      !!
70      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
71      !!----------------------------------------------------------------------
72      USE oce, ONLY:   zfu => ta,   & ! use ta as 3D workspace
73                       zfv => sa      ! use sa as 3D workspace
74
75      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
76
77      INTEGER  ::   ji, jj, jk                      ! dummy loop indices
78      REAL(wp) ::   zua, zva, zbu, zbv              ! temporary scalars
79      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v! temporary scalars
80      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f     ! temporary workspace
81      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f     !    "           "
82      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfw, zfu_uw, zfv_vw
83      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv   ! temporary workspace
84      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu   ! temporary workspace
85      !!----------------------------------------------------------------------
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      zfu_t(:,:,:) = 0.e0
93      zfv_t(:,:,:) = 0.e0
94      zfu_f(:,:,:) = 0.e0
95      zfv_f(:,:,:) = 0.e0
96     
97      zlu_uu(:,:,:,:) = 0.e0 
98      zlv_vv(:,:,:,:) = 0.e0 
99      zlu_uv(:,:,:,:) = 0.e0 
100      zlv_vu(:,:,:,:) = 0.e0 
101
102      IF( l_trddyn ) THEN           ! Save ua and va trends
103         zfu_uw(:,:,:) = ua(:,:,:)
104         zfv_vw(:,:,:) = va(:,:,:)
105      ENDIF
106
107      !                                                ! ===============
108      DO jk = 1, jpkm1                                 ! Horizontal slab
109         !                                             ! ===============
110         ! Laplacian
111         ! ---------
112         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
113         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
114           
115         DO jj = 2, jpjm1
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         !                                             ! ===============
129      END DO                                           !   End of slab
130      !                                                ! ===============
131
132      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.) 
133      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.) 
134      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.) 
135      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.) 
136
137      CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) 
138      CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) 
139      CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) 
140      CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 
141
142      ! I. Horizontal advection
143      ! -----------------------
144      !                                                ! ===============
145      DO jk = 1, jpkm1                                 ! Horizontal slab
146         !                                             ! ===============
147         ! horizontal volume fluxes
148         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
149         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
150
151         ! horizontal momentum fluxes at T- and F-point
152         DO jj = 1, jpjm1
153            DO ji = 1, fs_jpim1   ! vector opt.
154               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
155               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
156
157               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
158               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
159               ENDIF
160               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
161               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
162               ENDIF
163
164               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
165                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
166                  &                * ( zui - gamma1 * zl_u)
167               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
168                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
169                  &                * ( zvj - gamma1 * zl_v)
170
171               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
172               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
173               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
174               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
175               ENDIF
176               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
177               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
178               ENDIF
179
180               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
181                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
182               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
183                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
184            END DO
185         END DO
186
187         ! divergence of horizontal momentum fluxes
188         DO jj = 2, jpjm1
189            DO ji = fs_2, fs_jpim1   ! vector opt.
190               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
191               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
192               ! horizontal advective trends
193               zua = - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)   &
194                  &     + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
195               zva = - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)   &
196                  &     + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
197               ! add it to the general tracer trends
198               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
199               va(ji,jj,jk) = va(ji,jj,jk) + zva
200            END DO
201         END DO
202         !                                             ! ===============
203      END DO                                           !   End of slab
204      !                                                ! ===============
205
206      IF( l_trddyn ) THEN      ! save the horizontal advection trend for diagnostic
207         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
208         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
209         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
210      ENDIF
211
212      ! II. Vertical advection
213      ! ----------------------
214
215      IF( l_trddyn ) THEN           ! Save ua and va trends
216         zfu_t(:,:,:) = ua(:,:,:)
217         zfv_t(:,:,:) = va(:,:,:)
218      ENDIF
219
220      ! Second order centered tracer flux at w-point
221      DO jk = 1, jpkm1
222         ! Vertical volume fluxes                   
223         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
224         ! Vertical advective fluxes                   
225         IF( jk == 1 ) THEN
226            zfu_uw(:,:,jpk) = 0.e0         ! Bottom value : flux set to zero
227            zfv_vw(:,:,jpk) = 0.e0
228            !                              ! Surface value
229            IF( lk_dynspg_rl ) THEN             ! rigid lid : flux set to zero
230               zfu_uw(:,:, 1 ) = 0.e0   
231               zfv_vw(:,:, 1 ) = 0.e0
232            ELSE                                ! free surface-constant volume
233               DO jj = 2, jpjm1
234                  DO ji = fs_2, fs_jpim1
235                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
236                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
237                  END DO
238               END DO
239            ENDIF
240         ELSE
241            !                              ! interior fluxes
242            DO jj = 2, jpjm1
243               DO ji = fs_2, fs_jpim1   ! vector opt.
244                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
245                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
246               END DO
247            END DO
248         ENDIF
249      END DO
250
251      ! momentum flux divergence at u-, v-points added to the general trend
252      DO jk = 1, jpkm1
253         DO jj = 2, jpjm1 
254            DO ji = fs_2, fs_jpim1   ! vector opt.
255               zua = - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
256                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
257               zva = - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
258                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
259               ! add it to the general tracer trends
260               ua(ji,jj,jk) =  ua(ji,jj,jk) + zua
261               va(ji,jj,jk) =  va(ji,jj,jk) + zva
262            END DO
263         END DO
264      END DO
265
266      IF( l_trddyn ) THEN      ! save the vertical advection trend for diagnostic
267         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
268         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
269         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
270      ENDIF
271
272      !                             ! Control print
273      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
274         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
275
276   END SUBROUTINE dyn_adv_ubs
277
278   !!==============================================================================
279END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.