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/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN/dynadv_ubs.F90 @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago

Merge trunk -r14642:HEAD
  • Property svn:keywords set to Id
File size: 13.2 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"
[13237]36#  include "domzgr_substitute.h90"
[14219]37#  include "single_precision_substitute.h90"
[643]38   !!----------------------------------------------------------------------
[9598]39   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]40   !! $Id$
[10068]41   !! Software governed by the CeCILL license (see ./LICENSE)
[643]42   !!----------------------------------------------------------------------
43CONTAINS
44
[12377]45   SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs )
[643]46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE dyn_adv_ubs  ***
48      !!
49      !! ** Purpose :   Compute the now momentum advection trend in flux form
[1566]50      !!              and the general trend of the momentum equation.
[643]51      !!
52      !! ** Method  :   The scheme is the one implemeted in ROMS. It depends
53      !!      on two parameter gamma1 and gamma2. The former control the
54      !!      upstream baised part of the scheme and the later the centred
55      !!      part:     gamma1 = 0    pure centered  (no diffusive part)
56      !!                       = 1/4  Quick scheme
57      !!                       = 1/3  3rd order Upstream biased scheme
58      !!                gamma2 = 0    2nd order finite differencing
[4153]59      !!                       = 1/32 4th order finite differencing
[643]60      !!      For stability reasons, the first term of the fluxes which cor-
61      !!      responds to a second order centered scheme is evaluated using 
62      !!      the now velocity (centered in time) while the second term which 
63      !!      is the diffusive part of the scheme, is evaluated using the
64      !!      before velocity (forward in time).
65      !!      Default value (hard coded in the begining of the module) are
[4153]66      !!      gamma1=1/3 and gamma2=1/32.
[643]67      !!
[12377]68      !! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the 3D advective momentum trends
[643]69      !!
70      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
71      !!----------------------------------------------------------------------
[12377]72      INTEGER                             , INTENT( in )  ::  kt              ! ocean time-step index
73      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs  ! ocean time level indices
[14219]74      REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv        ! ocean velocities and RHS of momentum equation
[2715]75      !
[6140]76      INTEGER  ::   ji, jj, jk   ! dummy loop indices
77      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! local scalars
[14219]78      REAL(dp), DIMENSION(jpi,jpj,jpk)   :: zfu_t, zfu_uw
79      REAL(wp), DIMENSION(jpi,jpj,jpk)   :: zfu_f, zfu
80      REAL(dp), DIMENSION(jpi,jpj,jpk)    :: zfv_t, zfv_vw
81      REAL(wp), DIMENSION(jpi,jpj,jpk)   :: zfv_f, zfv, zfw
[9019]82      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv
83      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu
[643]84      !!----------------------------------------------------------------------
[3294]85      !
[643]86      IF( kt == nit000 ) THEN
87         IF(lwp) WRITE(numout,*)
88         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
89         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
90      ENDIF
[3294]91      !
[2715]92      zfu_t(:,:,:) = 0._wp
93      zfv_t(:,:,:) = 0._wp
94      zfu_f(:,:,:) = 0._wp
95      zfv_f(:,:,:) = 0._wp
[1566]96      !
[2715]97      zlu_uu(:,:,:,:) = 0._wp
98      zlv_vv(:,:,:,:) = 0._wp 
99      zlu_uv(:,:,:,:) = 0._wp 
100      zlv_vu(:,:,:,:) = 0._wp 
[6140]101      !
102      IF( l_trddyn ) THEN           ! trends: store the input trends
[12377]103         zfu_uw(:,:,:) = puu(:,:,:,Krhs)
104         zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
[1129]105      ENDIF
[1566]106      !                                      ! =========================== !
107      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
108         !                                   ! =========================== !
109         !                                         ! horizontal volume fluxes
[12377]110         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
111         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
[1566]112         !           
[13497]113         DO_2D( 0, 0, 0, 0 )                       ! laplacian
[12377]114            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)
115            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)
116            zlu_uv(ji,jj,jk,1) = ( puu (ji  ,jj+1,jk,Kbb) - puu (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
117               &               - ( puu (ji  ,jj  ,jk,Kbb) - puu (ji  ,jj-1,jk,Kbb) ) * fmask(ji  ,jj-1,jk)
118            zlv_vu(ji,jj,jk,1) = ( pvv (ji+1,jj  ,jk,Kbb) - pvv (ji  ,jj  ,jk,Kbb) ) * fmask(ji  ,jj  ,jk)   &
119               &               - ( pvv (ji  ,jj  ,jk,Kbb) - pvv (ji-1,jj  ,jk,Kbb) ) * fmask(ji-1,jj  ,jk)
120            !
121            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)
122            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)
123            zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
124               &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
125            zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
126               &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
127         END_2D
[1566]128      END DO
[14644]129      CALL lbc_lnk( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1.0_wp , zlu_uv(:,:,:,1), 'U', 1.0_wp,  &
130         &                        zlu_uu(:,:,:,2), 'U', 1.0_wp , zlu_uv(:,:,:,2), 'U', 1.0_wp,  & 
131         &                        zlv_vv(:,:,:,1), 'V', 1.0_wp , zlv_vu(:,:,:,1), 'V', 1.0_wp,  &
132         &                        zlv_vv(:,:,:,2), 'V', 1.0_wp , zlv_vu(:,:,:,2), 'V', 1.0_wp   )
[6140]133      !
[1566]134      !                                      ! ====================== !
135      !                                      !  Horizontal advection  !
136      DO jk = 1, jpkm1                       ! ====================== !
137         !                                         ! horizontal volume fluxes
[12377]138         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm)
139         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm)
[1566]140         !
[13497]141         DO_2D( 1, 0, 1, 0 )                       ! horizontal momentum fluxes at T- and F-point
[12377]142            zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj  ,jk,Kmm) )
143            zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji  ,jj+1,jk,Kmm) )
144            !
145            IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
146            ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
147            ENDIF
148            IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
149            ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
150            ENDIF
151            !
152            zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
153               &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
154               &                * ( zui - gamma1 * zl_u)
155            zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
156               &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
157               &                * ( zvj - gamma1 * zl_v)
158            !
159            zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
160            zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
161            IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
162            ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
163            ENDIF
164            IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
165            ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
166            ENDIF
167            !
168            zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
169               &                * ( puu(ji,jj,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) - gamma1 * zl_u )
170            zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
171               &                * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) - gamma1 * zl_v )
172         END_2D
[13497]173         DO_2D( 0, 0, 0, 0 )                       ! divergence of horizontal momentum fluxes
[12377]174            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
[13237]175               &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj)   &
176               &                           / e3u(ji,jj,jk,Kmm)
[12377]177            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
[13237]178               &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj)   &
179               &                           / e3v(ji,jj,jk,Kmm)
[12377]180         END_2D
[1566]181      END DO
[6140]182      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
[12377]183         zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
184         zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
185         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
186         zfu_t(:,:,:) = puu(:,:,:,Krhs)
187         zfv_t(:,:,:) = pvv(:,:,:,Krhs)
[1129]188      ENDIF
[1566]189      !                                      ! ==================== !
190      !                                      !  Vertical advection  !
[6140]191      !                                      ! ==================== !
[13497]192      DO_2D( 0, 0, 0, 0 )                          ! surface/bottom advective fluxes set to zero
[12377]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_2D
[6140]198      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
[13295]199         DO_2D( 0, 0, 0, 0 )
[12377]200            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)
201            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)
202         END_2D
[6140]203      ENDIF
204      DO jk = 2, jpkm1                          ! interior fluxes
[13295]205         DO_2D( 0, 1, 0, 1 )
[12377]206            zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
207         END_2D
[13295]208         DO_2D( 0, 0, 0, 0 )
[12377]209            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) )
210            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) )
211         END_2D
[643]212      END DO
[13497]213      DO_3D( 0, 0, 0, 0, 1, jpkm1 )             ! divergence of vertical momentum flux divergence
[13237]214         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)   &
215            &                                       / e3u(ji,jj,jk,Kmm)
216         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)   &
217            &                                       / e3v(ji,jj,jk,Kmm)
[12377]218      END_3D
[1566]219      !
[6140]220      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
[12377]221         zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
222         zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
223         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
[1129]224      ENDIF
[6140]225      !                                         ! Control print
[14219]226      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=CASTWP(puu(:,:,:,Krhs)), clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
227         &                                  tab3d_2=CASTWP(pvv(:,:,:,Krhs)), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]228      !
[643]229   END SUBROUTINE dyn_adv_ubs
230
231   !!==============================================================================
232END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.