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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 13321

Last change on this file since 13321 was 11101, checked in by frrh, 5 years ago

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