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 NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DYN/dynadv_ubs.F90

Last change on this file was 10986, checked in by andmirek, 5 years ago

GMED 462 add flush
File size: 12.9 KB
RevLine 
[643]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   !!======================================================================
[1566]7   !! History :  2.0  ! 2006-08  (R. Benshila, L. Debreu)  Original code
8   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
[643]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
[5062]18   USE trd_oce        ! trends: ocean variables
19   USE trddyn         ! trend manager: dynamics
20   !
[2715]21   USE in_out_manager ! I/O manager
[1129]22   USE prtctl         ! Print control
[2715]23   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
24   USE lib_mpp        ! MPP library
[643]25
26   IMPLICIT NONE
27   PRIVATE
28
[3294]29   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS
[4153]30   REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0   2nd order  ; =1/32 4th order centred
[643]31
[1566]32   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
[643]33
34   !! * Substitutions
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
[9598]37   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]38   !! $Id$
[10068]39   !! Software governed by the CeCILL license (see ./LICENSE)
[643]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
[1566]48      !!              and the general trend of the momentum equation.
[643]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
[4153]57      !!                       = 1/32 4th order finite differencing
[643]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
[4153]64      !!      gamma1=1/3 and gamma2=1/32.
[643]65      !!
[1566]66      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
[643]67      !!
68      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
69      !!----------------------------------------------------------------------
[6140]70      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
[2715]71      !
[6140]72      INTEGER  ::   ji, jj, jk   ! dummy loop indices
73      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! local scalars
[9019]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
[643]78      !!----------------------------------------------------------------------
[3294]79      !
[643]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,*) '~~~~~~~~~~~'
[10986]84         IF(lflush) CALL FLUSH(numout)
[643]85      ENDIF
[3294]86      !
[2715]87      zfu_t(:,:,:) = 0._wp
88      zfv_t(:,:,:) = 0._wp
89      zfu_f(:,:,:) = 0._wp
90      zfv_f(:,:,:) = 0._wp
[1566]91      !
[2715]92      zlu_uu(:,:,:,:) = 0._wp
93      zlv_vv(:,:,:,:) = 0._wp 
94      zlu_uv(:,:,:,:) = 0._wp 
95      zlv_vu(:,:,:,:) = 0._wp 
[6140]96      !
97      IF( l_trddyn ) THEN           ! trends: store the input trends
[1129]98         zfu_uw(:,:,:) = ua(:,:,:)
99         zfv_vw(:,:,:) = va(:,:,:)
100      ENDIF
[1566]101      !                                      ! =========================== !
102      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
103         !                                   ! =========================== !
104         !                                         ! horizontal volume fluxes
[6140]105         zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
106         zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
[1566]107         !           
108         DO jj = 2, jpjm1                          ! laplacian
[643]109            DO ji = fs_2, fs_jpim1   ! vector opt.
[5069]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) - ub (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
113                  &               - ( ub (ji  ,jj  ,jk) - ub (ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
114               zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj  ,jk) - vb (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
115                  &               - ( vb (ji  ,jj  ,jk) - vb (ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
[2715]116               !
[5069]117               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)
118               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)
119               zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
120                  &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
121               zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
122                  &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
[643]123            END DO
124         END DO
[1566]125      END DO
[10425]126      CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1. , zlu_uv(:,:,:,1), 'U', 1.,  &
[9090]127                      &   zlu_uu(:,:,:,2), 'U', 1. , zlu_uv(:,:,:,2), 'U', 1.,  & 
128                      &   zlv_vv(:,:,:,1), 'V', 1. , zlv_vu(:,:,:,1), 'V', 1.,  &
129                      &   zlv_vv(:,:,:,2), 'V', 1. , zlv_vu(:,:,:,2), 'V', 1.   )
[6140]130      !
[1566]131      !                                      ! ====================== !
132      !                                      !  Horizontal advection  !
133      DO jk = 1, jpkm1                       ! ====================== !
134         !                                         ! horizontal volume fluxes
[6140]135         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
136         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
[1566]137         !
138         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
[643]139            DO ji = 1, fs_jpim1   ! vector opt.
140               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
141               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
[1566]142               !
[6140]143               IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
144               ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
[643]145               ENDIF
[6140]146               IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
147               ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
[643]148               ENDIF
[1566]149               !
[643]150               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
151                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
152                  &                * ( zui - gamma1 * zl_u)
153               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
154                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
155                  &                * ( zvj - gamma1 * zl_v)
[1566]156               !
[643]157               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
158               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
[6140]159               IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
160               ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
[643]161               ENDIF
[6140]162               IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
163               ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
[643]164               ENDIF
[1566]165               !
[643]166               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
167                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
168               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
169                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
170            END DO
171         END DO
[1566]172         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
[643]173            DO ji = fs_2, fs_jpim1   ! vector opt.
[6140]174               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
175                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
176               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
177                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
[643]178            END DO
179         END DO
[1566]180      END DO
[6140]181      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
[1129]182         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
183         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
[5062]184         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
[1129]185         zfu_t(:,:,:) = ua(:,:,:)
186         zfv_t(:,:,:) = va(:,:,:)
187      ENDIF
[1566]188      !                                      ! ==================== !
189      !                                      !  Vertical advection  !
[6140]190      !                                      ! ==================== !
191      DO jj = 2, jpjm1                             ! surface/bottom advective fluxes set to zero                 
192         DO ji = fs_2, fs_jpim1
193            zfu_uw(ji,jj,jpk) = 0._wp
194            zfv_vw(ji,jj,jpk) = 0._wp
195            zfu_uw(ji,jj, 1 ) = 0._wp
196            zfv_vw(ji,jj, 1 ) = 0._wp
197         END DO
198      END DO
199      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
200         DO jj = 2, jpjm1
201            DO ji = fs_2, fs_jpim1
202               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)
203               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)
[643]204            END DO
[6140]205         END DO
206      ENDIF
207      DO jk = 2, jpkm1                          ! interior fluxes
[6750]208         DO jj = 2, jpj
209            DO ji = 2, jpi
[6140]210               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk)
211            END DO
212         END DO
213         DO jj = 2, jpjm1
214            DO ji = fs_2, fs_jpim1   ! vector opt.
215               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
216               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
217            END DO
218         END DO
[643]219      END DO
[6140]220      DO jk = 1, jpkm1                          ! divergence of vertical momentum flux divergence
221         DO jj = 2, jpjm1
[643]222            DO ji = fs_2, fs_jpim1   ! vector opt.
[6140]223               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)
224               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)
[643]225            END DO
226         END DO
227      END DO
[1566]228      !
[6140]229      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
[1129]230         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
231         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
[5062]232         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
[1129]233      ENDIF
[6140]234      !                                         ! Control print
[1129]235      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
236         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]237      !
[643]238   END SUBROUTINE dyn_adv_ubs
239
240   !!==============================================================================
241END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.