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/2019/fix_sn_cfctl_ticket2328/src/OCE/DYN – NEMO

source: NEMO/branches/2019/fix_sn_cfctl_ticket2328/src/OCE/DYN/dynadv_ubs.F90 @ 11869

Last change on this file since 11869 was 11869, checked in by acc, 4 years ago

Branch 2019/fix_sn_cfctl_ticket2328. Changes to enable correct functionality for the sn_cfctl%l_mppout and sn_cfctl%l_mpptop options. These changes also introduce a sn_cfctl%l_oasout option to toggle the OASIS setup information (sbccpl.F90, only) which was yet another misuse of ln_ctl. The next step may be to remove most references to ln_ctl altogether and provide a single control mechanism. TBD. See ticket #2328
  • Property svn:keywords set to Id
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,*) '~~~~~~~~~~~'
84      ENDIF
[3294]85      !
[2715]86      zfu_t(:,:,:) = 0._wp
87      zfv_t(:,:,:) = 0._wp
88      zfu_f(:,:,:) = 0._wp
89      zfv_f(:,:,:) = 0._wp
[1566]90      !
[2715]91      zlu_uu(:,:,:,:) = 0._wp
92      zlv_vv(:,:,:,:) = 0._wp 
93      zlu_uv(:,:,:,:) = 0._wp 
94      zlv_vu(:,:,:,:) = 0._wp 
[6140]95      !
96      IF( l_trddyn ) THEN           ! trends: store the input trends
[1129]97         zfu_uw(:,:,:) = ua(:,:,:)
98         zfv_vw(:,:,:) = va(:,:,:)
99      ENDIF
[1566]100      !                                      ! =========================== !
101      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
102         !                                   ! =========================== !
103         !                                         ! horizontal volume fluxes
[6140]104         zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
105         zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
[1566]106         !           
107         DO jj = 2, jpjm1                          ! laplacian
[643]108            DO ji = fs_2, fs_jpim1   ! vector opt.
[5069]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)
[2715]115               !
[5069]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)
[643]122            END DO
123         END DO
[1566]124      END DO
[10425]125      CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1. , zlu_uv(:,:,:,1), 'U', 1.,  &
[9090]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.   )
[6140]129      !
[1566]130      !                                      ! ====================== !
131      !                                      !  Horizontal advection  !
132      DO jk = 1, jpkm1                       ! ====================== !
133         !                                         ! horizontal volume fluxes
[6140]134         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk)
135         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk)
[1566]136         !
137         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
[643]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) )
[1566]141               !
[6140]142               IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
143               ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
[643]144               ENDIF
[6140]145               IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
146               ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
[643]147               ENDIF
[1566]148               !
[643]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)
[1566]155               !
[643]156               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
157               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
[6140]158               IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
159               ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
[643]160               ENDIF
[6140]161               IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
162               ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
[643]163               ENDIF
[1566]164               !
[643]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
[1566]171         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
[643]172            DO ji = fs_2, fs_jpim1   ! vector opt.
[6140]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)
[643]177            END DO
178         END DO
[1566]179      END DO
[6140]180      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
[1129]181         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
182         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
[5062]183         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
[1129]184         zfu_t(:,:,:) = ua(:,:,:)
185         zfv_t(:,:,:) = va(:,:,:)
186      ENDIF
[1566]187      !                                      ! ==================== !
188      !                                      !  Vertical advection  !
[6140]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)
[643]203            END DO
[6140]204         END DO
205      ENDIF
206      DO jk = 2, jpkm1                          ! interior fluxes
[6750]207         DO jj = 2, jpj
208            DO ji = 2, jpi
[6140]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
[643]218      END DO
[6140]219      DO jk = 1, jpkm1                          ! divergence of vertical momentum flux divergence
220         DO jj = 2, jpjm1
[643]221            DO ji = fs_2, fs_jpim1   ! vector opt.
[6140]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)
[643]224            END DO
225         END DO
226      END DO
[1566]227      !
[6140]228      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
[1129]229         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
230         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
[5062]231         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
[1129]232      ENDIF
[6140]233      !                                         ! Control print
[11869]234      IF(ln_ctl .OR. sn_cfctl%l_mppout)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
235         &                                              tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]236      !
[643]237   END SUBROUTINE dyn_adv_ubs
238
239   !!==============================================================================
240END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.