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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 2620

Last change on this file since 2620 was 2590, checked in by trackstand2, 13 years ago

Merge branch 'dynamic_memory' into master-svn-dyn
  • Property svn:keywords set to Id
File size: 14.9 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 in_out_manager ! I/O manager
19   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
20   USE trdmod         ! ocean dynamics trends
21   USE trdmod_oce     ! ocean variables trends
22   USE prtctl         ! Print control
23
24   IMPLICIT NONE
25   PRIVATE
26
27   REAL(wp), PARAMETER :: gamma1 = 1._wp/4._wp  ! =1/4 quick      ; =1/3  3rd order UBS
28   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred
29
30   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
37   !! $Id$
38   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
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 oce, ONLY:   zfu => ta   ! use ta as 3D workspace
71      USE oce, ONLY:   zfv => sa   ! use sa as 3D workspace
72      USE wrk_nemo, ONLY: wrk_use, wrk_release
73      USE wrk_nemo, ONLY: zfu_t  =>wrk_3d_1, &
74                          zfu_f  =>wrk_3d_2, &
75                          zfv_t  =>wrk_3d_3, &
76                          zfv_f  =>wrk_3d_4, &
77                          zfw    =>wrk_3d_5, &
78                          zfu_uw =>wrk_3d_6, &
79                          zfv_vw =>wrk_3d_7 
80      USE wrk_nemo, ONLY: zlu_uu=>wrk_4d_1, &
81                          zlu_uv=>wrk_4d_2, &
82                          zlv_vv=>wrk_4d_3, &
83                          zlv_vu=>wrk_4d_4
84      !!
85      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
86      !!
87      INTEGER  ::   ji, jj, jk            ! dummy loop indices
88      REAL(wp) ::   zbu, zbv    ! temporary scalars
89      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! temporary scalars
90! ARPDBG - arrays below replaced with global work spaces
91!!$      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f     ! temporary workspace
92!!$      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f     !    "           "
93!!$      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfw, zfu_uw, zfv_vw
94!!$      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv   ! temporary workspace
95!!$      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu   ! temporary workspace
96      !!----------------------------------------------------------------------
97
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
101         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
102      ENDIF
103
104      ! Check that required workspace arrays are not already in use
105      IF( .not. wrk_use(3, 1, 2, 3, 4, 5, 6, 7) )THEN
106         CALL ctl_stop('dyn_adv_ubs : error : required 3d workspace array is already in use')
107      END IF
108      IF(.not. wrk_use(4, 1, 2, 3, 4) )THEN
109         CALL ctl_stop('dyn_adv_ubs : error : required 4d workspace array is already in use')
110      END IF
111
112      zfu_t(:,:,:) = 0.e0
113      zfv_t(:,:,:) = 0.e0
114      zfu_f(:,:,:) = 0.e0
115      zfv_f(:,:,:) = 0.e0
116      !
117      zlu_uu(:,:,:,:) = 0.e0 
118      zlv_vv(:,:,:,:) = 0.e0 
119      zlu_uv(:,:,:,:) = 0.e0 
120      zlv_vu(:,:,:,:) = 0.e0 
121
122      IF( l_trddyn ) THEN           ! Save ua and va trends
123         zfu_uw(:,:,:) = ua(:,:,:)
124         zfv_vw(:,:,:) = va(:,:,:)
125      ENDIF
126
127      !                                      ! =========================== !
128      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
129         !                                   ! =========================== !
130         !                                         ! horizontal volume fluxes
131         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
132         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
133         !           
134         DO jj = 2, jpjm1                          ! laplacian
135            DO ji = fs_2, fs_jpim1   ! vector opt.
136               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)
137               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) 
138               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)
139               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)
140               
141               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)
142               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)
143               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)
144               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)
145            END DO
146         END DO
147      END DO
148!!gm BUG !!!  just below this should be +1 in all the communications
149      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.)
150      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.)
151      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.)
152      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 
153
154!!gm corrected:
155      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
156      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
157      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
158      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 
159!!gm end
160     
161      !                                      ! ====================== !
162      !                                      !  Horizontal advection  !
163      DO jk = 1, jpkm1                       ! ====================== !
164         !                                         ! horizontal volume fluxes
165         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
166         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
167         !
168         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
169            DO ji = 1, fs_jpim1   ! vector opt.
170               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
171               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
172               !
173               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
174               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
175               ENDIF
176               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
177               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
178               ENDIF
179               !
180               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
181                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
182                  &                * ( zui - gamma1 * zl_u)
183               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
184                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
185                  &                * ( zvj - gamma1 * zl_v)
186               !
187               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
188               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
189               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
190               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
191               ENDIF
192               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
193               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
194               ENDIF
195               !
196               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
197                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
198               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
199                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
200            END DO
201         END DO
202         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
203            DO ji = fs_2, fs_jpim1   ! vector opt.
204               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
205               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
206               !
207               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
208                  &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
209               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    &
210                  &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
211            END DO
212         END DO
213      END DO
214      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
215         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
216         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
217         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
218         zfu_t(:,:,:) = ua(:,:,:)
219         zfv_t(:,:,:) = va(:,:,:)
220      ENDIF
221
222      !                                      ! ==================== !
223      !                                      !  Vertical advection  !
224      DO jk = 1, jpkm1                       ! ==================== !
225         !                                         ! Vertical volume fluxesÊ
226         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
227         !
228         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
229            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
230            zfv_vw(:,:,jpk) = 0.e0
231            !                                           ! Surface value :
232            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
233               zfu_uw(:,:, 1 ) = 0.e0   
234               zfv_vw(:,:, 1 ) = 0.e0
235            ELSE                                             ! constant volume : advection through the surface
236               DO jj = 2, jpjm1
237                  DO ji = fs_2, fs_jpim1
238                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
239                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
240                  END DO
241               END DO
242            ENDIF
243         ELSE                                      ! interior fluxes
244            DO jj = 2, jpjm1
245               DO ji = fs_2, fs_jpim1   ! vector opt.
246                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
247                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
248               END DO
249            END DO
250         ENDIF
251      END DO
252      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
253         DO jj = 2, jpjm1 
254            DO ji = fs_2, fs_jpim1   ! vector opt.
255               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
256                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
257               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
258                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
259            END DO
260         END DO
261      END DO
262      !
263      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
264         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
265         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
266         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
267      ENDIF
268      !                                            ! Control print
269      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
270         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
271      !
272      ! Signal that we're done with the 3D and 4D global workspace arrays
273      IF( (.not. wrk_release(3, 1, 2, 3, 4, 5, 6, 7)) .OR. &
274          (.not. wrk_release(4, 1, 2, 3, 4)) )THEN
275         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : failed to release workspace arrays'
276      END IF
277      !
278   END SUBROUTINE dyn_adv_ubs
279
280   !!==============================================================================
281END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.