source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90 @ 10874

Last change on this file since 10874 was 10874, checked in by davestorkey, 20 months ago

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Revert all changes so far in preparation for implementation of new design.

  • Property svn:keywords set to Id
File size: 12.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 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
26   IMPLICIT NONE
27   PRIVATE
28
29   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS
30   REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0   2nd order  ; =1/32 4th order centred
31
32   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
33
34   !! * Substitutions
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
38   !! $Id$
39   !! Software governed by the CeCILL license (see ./LICENSE)
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/32 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/3 and gamma2=1/32.
65      !!
66      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
67      !!
68      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
69      !!----------------------------------------------------------------------
70      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
71      !
72      INTEGER  ::   ji, jj, jk   ! dummy loop indices
73      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! local scalars
74      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f, zfu_uw, zfu
75      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f, zfv_vw, zfv, zfw
76      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv
77      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu
78      !!----------------------------------------------------------------------
79      !
80      IF( kt == nit000 ) THEN
81         IF(lwp) WRITE(numout,*)
82         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
83         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
84      ENDIF
85      !
86      zfu_t(:,:,:) = 0._wp
87      zfv_t(:,:,:) = 0._wp
88      zfu_f(:,:,:) = 0._wp
89      zfv_f(:,:,:) = 0._wp
90      !
91      zlu_uu(:,:,:,:) = 0._wp
92      zlv_vv(:,:,:,:) = 0._wp 
93      zlu_uv(:,:,:,:) = 0._wp 
94      zlv_vu(:,:,:,:) = 0._wp 
95      !
96      IF( l_trddyn ) THEN           ! trends: store the input trends
97         zfu_uw(:,:,:) = ua(:,:,:)
98         zfv_vw(:,:,:) = va(:,:,:)
99      ENDIF
100      !                                      ! =========================== !
101      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
102         !                                   ! =========================== !
103         !                                         ! horizontal volume fluxes
104         zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
105         zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
106         !           
107         DO jj = 2, jpjm1                          ! laplacian
108            DO ji = fs_2, fs_jpim1   ! vector opt.
109               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)
110               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)
111               zlu_uv(ji,jj,jk,1) = ( ub (ji  ,jj+1,jk) - ub (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
112                  &               - ( ub (ji  ,jj  ,jk) - ub (ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
113               zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj  ,jk) - vb (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
114                  &               - ( vb (ji  ,jj  ,jk) - vb (ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
115               !
116               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)
117               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)
118               zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
119                  &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
120               zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
121                  &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
122            END DO
123         END DO
124      END DO
125      CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1. , zlu_uv(:,:,:,1), 'U', 1.,  &
126                      &   zlu_uu(:,:,:,2), 'U', 1. , zlu_uv(:,:,:,2), 'U', 1.,  & 
127                      &   zlv_vv(:,:,:,1), 'V', 1. , zlv_vu(:,:,:,1), 'V', 1.,  &
128                      &   zlv_vv(:,:,:,2), 'V', 1. , zlv_vu(:,:,:,2), 'V', 1.   )
129      !
130      !                                      ! ====================== !
131      !                                      !  Horizontal advection  !
132      DO jk = 1, jpkm1                       ! ====================== !
133         !                                         ! horizontal volume fluxes
134         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
135         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
136         !
137         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
138            DO ji = 1, fs_jpim1   ! vector opt.
139               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
140               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
141               !
142               IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
143               ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
144               ENDIF
145               IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
146               ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
147               ENDIF
148               !
149               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
150                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
151                  &                * ( zui - gamma1 * zl_u)
152               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
153                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
154                  &                * ( zvj - gamma1 * zl_v)
155               !
156               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
157               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
158               IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
159               ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
160               ENDIF
161               IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
162               ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
163               ENDIF
164               !
165               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
166                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
167               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
168                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
169            END DO
170         END DO
171         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
172            DO ji = fs_2, fs_jpim1   ! vector opt.
173               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
174                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
175               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
176                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
177            END DO
178         END DO
179      END DO
180      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
181         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
182         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
183         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
184         zfu_t(:,:,:) = ua(:,:,:)
185         zfv_t(:,:,:) = va(:,:,:)
186      ENDIF
187      !                                      ! ==================== !
188      !                                      !  Vertical advection  !
189      !                                      ! ==================== !
190      DO jj = 2, jpjm1                             ! surface/bottom advective fluxes set to zero                 
191         DO ji = fs_2, fs_jpim1
192            zfu_uw(ji,jj,jpk) = 0._wp
193            zfv_vw(ji,jj,jpk) = 0._wp
194            zfu_uw(ji,jj, 1 ) = 0._wp
195            zfv_vw(ji,jj, 1 ) = 0._wp
196         END DO
197      END DO
198      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
199         DO jj = 2, jpjm1
200            DO ji = fs_2, fs_jpim1
201               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1)
202               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1)
203            END DO
204         END DO
205      ENDIF
206      DO jk = 2, jpkm1                          ! interior fluxes
207         DO jj = 2, jpj
208            DO ji = 2, jpi
209               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
210            END DO
211         END DO
212         DO jj = 2, jpjm1
213            DO ji = fs_2, fs_jpim1   ! vector opt.
214               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
215               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
216            END DO
217         END DO
218      END DO
219      DO jk = 1, jpkm1                          ! divergence of vertical momentum flux divergence
220         DO jj = 2, jpjm1
221            DO ji = fs_2, fs_jpim1   ! vector opt.
222               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
223               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
224            END DO
225         END DO
226      END DO
227      !
228      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
229         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
230         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
231         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
232      ENDIF
233      !                                         ! Control print
234      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
235         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
236      !
237   END SUBROUTINE dyn_adv_ubs
238
239   !!==============================================================================
240END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.