source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 12 months ago

The Dr Hook changes from my perl code.

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