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

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 2789

Last change on this file since 2789 was 2789, checked in by cetlod, 13 years ago

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