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

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

get back to the nemo_v2_3 version for trunk

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 12.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      ! I/O manager
20   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
21
22   IMPLICIT NONE
23   PRIVATE
24
25   REAL(wp), PARAMETER :: gamma1 = 1._wp/4._wp  ! =1/4 quick      ; =1/3  3rd order UBS
26   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred
27
28   !! * Routine accessibility
29   PUBLIC dyn_adv_ubs                            ! routine called by step.F90
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!   OPA 9.0 , LODYC-IPSL  (2006)
36   !! $Header$
37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE dyn_adv_ubs( kt )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE dyn_adv_ubs  ***
45      !!
46      !! ** Purpose :   Compute the now momentum advection trend in flux form
47      !!      and the general trend of the momentum equation.
48      !!
49      !! ** Method  :   The scheme is the one implemeted in ROMS. It depends
50      !!      on two parameter gamma1 and gamma2. The former control the
51      !!      upstream baised part of the scheme and the later the centred
52      !!      part:     gamma1 = 0    pure centered  (no diffusive part)
53      !!                       = 1/4  Quick scheme
54      !!                       = 1/3  3rd order Upstream biased scheme
55      !!                gamma2 = 0    2nd order finite differencing
56      !!                       = 1/8  4th order finite differencing
57      !!      For stability reasons, the first term of the fluxes which cor-
58      !!      responds to a second order centered scheme is evaluated using 
59      !!      the now velocity (centered in time) while the second term which 
60      !!      is the diffusive part of the scheme, is evaluated using the
61      !!      before velocity (forward in time).
62      !!      Default value (hard coded in the begining of the module) are
63      !!      gamma1=1/4 and gamma2=1/8.
64      !!
65      !! ** Action : - update (ua,va) with the 3D advective momentum trends
66      !!
67      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
68      !!----------------------------------------------------------------------
69      USE oce, ONLY:   zfu => ta,   & ! use ta as 3D workspace
70                       zfv => sa      ! use sa as 3D workspace
71
72      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
73
74      INTEGER  ::   ji, jj, jk                      ! dummy loop indices
75      REAL(wp) ::   zua, zva, zbu, zbv              ! temporary scalars
76      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v! temporary scalars
77      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f     ! temporary workspace
78      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f     !    "           "
79      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfw, zfu_uw, zfv_vw
80      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv   ! temporary workspace
81      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu   ! temporary workspace
82      !!----------------------------------------------------------------------
83
84      IF( kt == nit000 ) THEN
85         IF(lwp) WRITE(numout,*)
86         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
87         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
88      ENDIF
89      zfu_t(:,:,:) = 0.e0
90      zfv_t(:,:,:) = 0.e0
91      zfu_f(:,:,:) = 0.e0
92      zfv_f(:,:,:) = 0.e0
93     
94      zlu_uu(:,:,:,:) = 0.e0 
95      zlv_vv(:,:,:,:) = 0.e0 
96      zlu_uv(:,:,:,:) = 0.e0 
97      zlv_vu(:,:,:,:) = 0.e0 
98
99     
100      !                                                ! ===============
101      DO jk = 1, jpkm1                                 ! Horizontal slab
102         !                                             ! ===============
103         ! Laplacian
104         ! ---------
105         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
106         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
107           
108         DO jj = 2, jpjm1
109            DO ji = fs_2, fs_jpim1   ! vector opt.
110               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)
111               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) 
112               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)
113               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)
114               
115               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)
116               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)
117               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)
118               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)
119            END DO
120         END DO
121         !                                             ! ===============
122      END DO                                           !   End of slab
123      !                                                ! ===============
124
125      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.) 
126      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.) 
127      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.) 
128      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.) 
129
130      CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) 
131      CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) 
132      CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) 
133      CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 
134
135      ! I. Horizontal advection
136      ! -----------------------
137      !                                                ! ===============
138      DO jk = 1, jpkm1                                 ! Horizontal slab
139         !                                             ! ===============
140         ! horizontal volume fluxes
141         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
142         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
143
144         ! horizontal momentum fluxes at T- and F-point
145         DO jj = 1, jpjm1
146            DO ji = 1, fs_jpim1   ! vector opt.
147               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
148               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
149
150               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
151               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
152               ENDIF
153               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
154               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
155               ENDIF
156
157               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
158                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
159                  &                * ( zui - gamma1 * zl_u)
160               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
161                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
162                  &                * ( zvj - gamma1 * zl_v)
163
164               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
165               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
166               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
167               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
168               ENDIF
169               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
170               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
171               ENDIF
172
173               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
174                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
175               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
176                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
177            END DO
178         END DO
179
180         ! divergence of horizontal momentum fluxes
181         DO jj = 2, jpjm1
182            DO ji = fs_2, fs_jpim1   ! vector opt.
183               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
184               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
185               ! horizontal advective trends
186               zua = - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)   &
187                  &     + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
188               zva = - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)   &
189                  &     + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
190               ! add it to the general tracer trends
191               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
192               va(ji,jj,jk) = va(ji,jj,jk) + zva
193#if defined key_trddyn
194               utrd(ji,jj,jk,1) = zua                    ! save the horizontal advective trend of momentum
195               vtrd(ji,jj,jk,1) = zva 
196#endif
197            END DO
198         END DO
199         !                                             ! ===============
200      END DO                                           !   End of slab
201      !                                                ! ===============
202
203
204      ! II. Vertical advection
205      ! ----------------------
206
207      ! Second order centered tracer flux at w-point
208      DO jk = 1, jpkm1
209         ! Vertical volume fluxes                   
210         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
211         ! Vertical advective fluxes                   
212         IF( jk == 1 ) THEN
213            zfu_uw(:,:,jpk) = 0.e0         ! Bottom value : flux set to zero
214            zfv_vw(:,:,jpk) = 0.e0
215            !                              ! Surface value
216            IF( lk_dynspg_rl ) THEN             ! rigid lid : flux set to zero
217               zfu_uw(:,:, 1 ) = 0.e0   
218               zfv_vw(:,:, 1 ) = 0.e0
219            ELSE                                ! free surface-constant volume
220               DO jj = 2, jpjm1
221                  DO ji = fs_2, fs_jpim1
222                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
223                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
224                  END DO
225               END DO
226            ENDIF
227         ELSE
228            !                              ! interior fluxes
229            DO jj = 2, jpjm1
230               DO ji = fs_2, fs_jpim1   ! vector opt.
231                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
232                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
233               END DO
234            END DO
235         ENDIF
236      END DO
237
238      ! momentum flux divergence at u-, v-points added to the general trend
239      DO jk = 1, jpkm1
240         DO jj = 2, jpjm1 
241            DO ji = fs_2, fs_jpim1   ! vector opt.
242               zua = - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
243                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
244               zva = - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
245                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
246               ! add it to the general tracer trends
247               ua(ji,jj,jk) =  ua(ji,jj,jk) + zua
248               va(ji,jj,jk) =  va(ji,jj,jk) + zva
249            END DO
250         END DO
251      END DO
252
253   END SUBROUTINE dyn_adv_ubs
254
255   !!==============================================================================
256END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.