source: NEMO/trunk/src/OCE/DYN/dynadv_ubs.F90 @ 13226

Last change on this file since 13226 was 13226, checked in by orioltp, 4 months ago

Merging dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation into the trunk

  • Property svn:keywords set to Id
File size: 12.5 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
[12377]35#  include "do_loop_substitute.h90"
[643]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
[12377]43   SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs )
[643]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      !!
[12377]66      !! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the 3D advective momentum trends
[643]67      !!
68      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
69      !!----------------------------------------------------------------------
[12377]70      INTEGER                             , INTENT( in )  ::  kt              ! ocean time-step index
71      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs  ! ocean time level indices
72      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv        ! ocean velocities and RHS of momentum equation
[2715]73      !
[6140]74      INTEGER  ::   ji, jj, jk   ! dummy loop indices
75      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! local scalars
[9019]76      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f, zfu_uw, zfu
77      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f, zfv_vw, zfv, zfw
78      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv
79      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu
[643]80      !!----------------------------------------------------------------------
[3294]81      !
[643]82      IF( kt == nit000 ) THEN
83         IF(lwp) WRITE(numout,*)
84         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
85         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
86      ENDIF
[3294]87      !
[2715]88      zfu_t(:,:,:) = 0._wp
89      zfv_t(:,:,:) = 0._wp
90      zfu_f(:,:,:) = 0._wp
91      zfv_f(:,:,:) = 0._wp
[1566]92      !
[2715]93      zlu_uu(:,:,:,:) = 0._wp
94      zlv_vv(:,:,:,:) = 0._wp 
95      zlu_uv(:,:,:,:) = 0._wp 
96      zlv_vu(:,:,:,:) = 0._wp 
[6140]97      !
98      IF( l_trddyn ) THEN           ! trends: store the input trends
[12377]99         zfu_uw(:,:,:) = puu(:,:,:,Krhs)
100         zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
[1129]101      ENDIF
[1566]102      !                                      ! =========================== !
103      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
104         !                                   ! =========================== !
105         !                                         ! horizontal volume fluxes
[12377]106         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
107         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
[1566]108         !           
[12377]109         DO_2D_00_00
110            zlu_uu(ji,jj,jk,1) = ( puu (ji+1,jj  ,jk,Kbb) - 2.*puu (ji,jj,jk,Kbb) + puu (ji-1,jj  ,jk,Kbb) ) * umask(ji,jj,jk)
111            zlv_vv(ji,jj,jk,1) = ( pvv (ji  ,jj+1,jk,Kbb) - 2.*pvv (ji,jj,jk,Kbb) + pvv (ji  ,jj-1,jk,Kbb) ) * vmask(ji,jj,jk)
112            zlu_uv(ji,jj,jk,1) = ( puu (ji  ,jj+1,jk,Kbb) - puu (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
113               &               - ( puu (ji  ,jj  ,jk,Kbb) - puu (ji  ,jj-1,jk,Kbb) ) * fmask(ji  ,jj-1,jk)
114            zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj  ,jk,Kbb) - pvv (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
115               &               - ( pvv (ji  ,jj  ,jk,Kbb) - pvv (ji-1,jj  ,jk,Kbb) ) * fmask(ji-1,jj  ,jk)
116            !
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)
123         END_2D
[1566]124      END DO
[13226]125      CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1.0_wp , zlu_uv(:,:,:,1), 'U', 1.0_wp,  &
126                      &   zlu_uu(:,:,:,2), 'U', 1.0_wp , zlu_uv(:,:,:,2), 'U', 1.0_wp,  & 
127                      &   zlv_vv(:,:,:,1), 'V', 1.0_wp , zlv_vu(:,:,:,1), 'V', 1.0_wp,  &
128                      &   zlv_vv(:,:,:,2), 'V', 1.0_wp , zlv_vu(:,:,:,2), 'V', 1.0_wp   )
[6140]129      !
[1566]130      !                                      ! ====================== !
131      !                                      !  Horizontal advection  !
132      DO jk = 1, jpkm1                       ! ====================== !
133         !                                         ! horizontal volume fluxes
[12377]134         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
135         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
[1566]136         !
[12377]137         DO_2D_10_10
138            zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
139            zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
140            !
141            IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
142            ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
143            ENDIF
144            IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
145            ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
146            ENDIF
147            !
148            zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
149               &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
150               &                * ( zui - gamma1 * zl_u)
151            zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
152               &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
153               &                * ( zvj - gamma1 * zl_v)
154            !
155            zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
156            zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
157            IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
158            ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
159            ENDIF
160            IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
161            ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
162            ENDIF
163            !
164            zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
165               &                * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) - gamma1 * zl_u )
166            zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
167               &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
168         END_2D
169         DO_2D_00_00
170            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
171               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
172            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
173               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
174         END_2D
[1566]175      END DO
[6140]176      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
[12377]177         zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
178         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
179         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
180         zfu_t(:,:,:) = puu(:,:,:,Krhs)
181         zfv_t(:,:,:) = pvv(:,:,:,Krhs)
[1129]182      ENDIF
[1566]183      !                                      ! ==================== !
184      !                                      !  Vertical advection  !
[6140]185      !                                      ! ==================== !
[12377]186      DO_2D_00_00
187         zfu_uw(ji,jj,jpk) = 0._wp
188         zfv_vw(ji,jj,jpk) = 0._wp
189         zfu_uw(ji,jj, 1 ) = 0._wp
190         zfv_vw(ji,jj, 1 ) = 0._wp
191      END_2D
[6140]192      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
[12377]193         DO_2D_00_00
194            zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
195            zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
196         END_2D
[6140]197      ENDIF
198      DO jk = 2, jpkm1                          ! interior fluxes
[12377]199         DO_2D_01_01
200            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
201         END_2D
202         DO_2D_00_00
203            zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
204            zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
205         END_2D
[643]206      END DO
[12377]207      DO_3D_00_00( 1, jpkm1 )
208         puu(ji,jj,jk,Krhs) =  puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
209         pvv(ji,jj,jk,Krhs) =  pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
210      END_3D
[1566]211      !
[6140]212      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
[12377]213         zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
214         zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
215         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
[1129]216      ENDIF
[6140]217      !                                         ! Control print
[12377]218      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
219         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]220      !
[643]221   END SUBROUTINE dyn_adv_ubs
222
223   !!==============================================================================
224END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.