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/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90 @ 10802

Last change on this file since 10802 was 10802, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : introduce new T/S variables and convert tracer advection routines (including calls from TOP).
  • Property svn:keywords set to Id
File size: 13.4 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
[10789]43   SUBROUTINE dyn_adv_ubs( kt, ktlev1, ktlev2, pu_rhs, pv_rhs )
[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      !!
[10789]66      !! ** Action : - (pu_rhs,pv_rhs) updated with the 3D advective momentum trends
[643]67      !!
68      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
69      !!----------------------------------------------------------------------
[10789]70      INTEGER, INTENT(in) ::   kt               ! ocean time-step index
71      INTEGER, INTENT(in) ::   ktlev1, ktlev2   ! time level indices for source terms
72      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends
[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
[10789]99         zfu_uw(:,:,:) = pu_rhs(:,:,:)
100         zfv_vw(:,:,:) = pv_rhs(:,:,:)
[1129]101      ENDIF
[1566]102      !                                      ! =========================== !
103      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
104         !                                   ! =========================== !
105         !                                         ! horizontal volume fluxes
[10789]106         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,ktlev2) * uu(:,:,jk,ktlev2)
107         zfv(:,:,jk) = e1v(:,:) * e3v(:,:,jk,ktlev2) * vv(:,:,jk,ktlev2)
[1566]108         !           
109         DO jj = 2, jpjm1                          ! laplacian
[643]110            DO ji = fs_2, fs_jpim1   ! vector opt.
[10789]111               zlu_uu(ji,jj,jk,1) = ( uu (ji+1,jj  ,jk,ktlev1) - 2.*uu (ji,jj,jk,ktlev1) + uu (ji-1,jj  ,jk,ktlev1) ) * umask(ji,jj,jk)
112               zlv_vv(ji,jj,jk,1) = ( vv (ji  ,jj+1,jk,ktlev1) - 2.*vv (ji,jj,jk,ktlev1) + vv (ji  ,jj-1,jk,ktlev1) ) * vmask(ji,jj,jk)
113               zlu_uv(ji,jj,jk,1) = ( uu (ji  ,jj+1,jk,ktlev1) - uu (ji  ,jj  ,jk,ktlev1) ) * fmask(ji  ,jj  ,jk)   &
114                  &               - ( uu (ji  ,jj  ,jk,ktlev1) - uu (ji  ,jj-1,jk,ktlev1) ) * fmask(ji  ,jj-1,jk)
115               zlv_vu(ji,jj,jk,1) = ( vv (ji+1,jj  ,jk,ktlev1) - vv (ji  ,jj  ,jk,ktlev1) ) * fmask(ji  ,jj  ,jk)   &
116                  &               - ( vv (ji  ,jj  ,jk,ktlev1) - vv (ji-1,jj  ,jk,ktlev1) ) * fmask(ji-1,jj  ,jk)
[2715]117               !
[5069]118               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)
119               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)
120               zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
121                  &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk)
122               zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   &
123                  &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk)
[643]124            END DO
125         END DO
[1566]126      END DO
[10425]127      CALL lbc_lnk_multi( 'dynadv_ubs', zlu_uu(:,:,:,1), 'U', 1. , zlu_uv(:,:,:,1), 'U', 1.,  &
[9090]128                      &   zlu_uu(:,:,:,2), 'U', 1. , zlu_uv(:,:,:,2), 'U', 1.,  & 
129                      &   zlv_vv(:,:,:,1), 'V', 1. , zlv_vu(:,:,:,1), 'V', 1.,  &
130                      &   zlv_vv(:,:,:,2), 'V', 1. , zlv_vu(:,:,:,2), 'V', 1.   )
[6140]131      !
[1566]132      !                                      ! ====================== !
133      !                                      !  Horizontal advection  !
134      DO jk = 1, jpkm1                       ! ====================== !
135         !                                         ! horizontal volume fluxes
[10789]136         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,ktlev2) * uu(:,:,jk,ktlev2)
137         zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v(:,:,jk,ktlev2) * vv(:,:,jk,ktlev2)
[1566]138         !
139         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
[643]140            DO ji = 1, fs_jpim1   ! vector opt.
[10789]141               zui = ( uu(ji,jj,jk,ktlev2) + uu(ji+1,jj  ,jk,ktlev2) )
142               zvj = ( vv(ji,jj,jk,ktlev2) + vv(ji  ,jj+1,jk,ktlev2) )
[1566]143               !
[6140]144               IF( zui > 0 ) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
145               ELSE                 ;   zl_u = zlu_uu(ji+1,jj,jk,1)
[643]146               ENDIF
[6140]147               IF( zvj > 0 ) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
148               ELSE                 ;   zl_v = zlv_vv(ji,jj+1,jk,1)
[643]149               ENDIF
[1566]150               !
[643]151               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
152                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
153                  &                * ( zui - gamma1 * zl_u)
154               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
155                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
156                  &                * ( zvj - gamma1 * zl_v)
[1566]157               !
[643]158               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
159               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
[6140]160               IF( zfuj > 0 ) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
161               ELSE                  ;    zl_v = zlv_vu( ji+1,jj,jk,1)
[643]162               ENDIF
[6140]163               IF( zfvi > 0 ) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
164               ELSE                  ;    zl_u = zlu_uv( ji,jj+1,jk,1)
[643]165               ENDIF
[1566]166               !
[643]167               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
[10789]168                  &                * ( uu(ji,jj,jk,ktlev2) + uu(ji  ,jj+1,jk,ktlev2) - gamma1 * zl_u )
[643]169               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
[10789]170                  &                * ( vv(ji,jj,jk,ktlev2) + vv(ji+1,jj  ,jk,ktlev2) - gamma1 * zl_v )
[643]171            END DO
172         END DO
[1566]173         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
[643]174            DO ji = fs_2, fs_jpim1   ! vector opt.
[10789]175               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - (  zfu_t(ji+1,jj,jk) - zfu_t(ji,jj  ,jk)    &
176                  &                           + zfv_f(ji  ,jj,jk) - zfv_f(ji,jj-1,jk)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev2)
177               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - (  zfu_f(ji,jj  ,jk) - zfu_f(ji-1,jj,jk)    &
178                  &                           + zfv_t(ji,jj+1,jk) - zfv_t(ji  ,jj,jk)  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev2)
[643]179            END DO
180         END DO
[1566]181      END DO
[6140]182      IF( l_trddyn ) THEN                          ! trends: send trends to trddyn for diagnostic
[10789]183         zfu_uw(:,:,:) = pu_rhs(:,:,:) - zfu_uw(:,:,:)
184         zfv_vw(:,:,:) = pv_rhs(:,:,:) - zfv_vw(:,:,:)
[5062]185         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
[10789]186         zfu_t(:,:,:) = pu_rhs(:,:,:)
187         zfv_t(:,:,:) = pv_rhs(:,:,:)
[1129]188      ENDIF
[1566]189      !                                      ! ==================== !
190      !                                      !  Vertical advection  !
[6140]191      !                                      ! ==================== !
192      DO jj = 2, jpjm1                             ! surface/bottom advective fluxes set to zero                 
193         DO ji = fs_2, fs_jpim1
194            zfu_uw(ji,jj,jpk) = 0._wp
195            zfv_vw(ji,jj,jpk) = 0._wp
196            zfu_uw(ji,jj, 1 ) = 0._wp
197            zfv_vw(ji,jj, 1 ) = 0._wp
198         END DO
199      END DO
200      IF( ln_linssh ) THEN                         ! constant volume : advection through the surface
201         DO jj = 2, jpjm1
202            DO ji = fs_2, fs_jpim1
[10789]203               zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * uu(ji,jj,1,ktlev2)
204               zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * vv(ji,jj,1,ktlev2)
[643]205            END DO
[6140]206         END DO
207      ENDIF
208      DO jk = 2, jpkm1                          ! interior fluxes
[6750]209         DO jj = 2, jpj
210            DO ji = 2, jpi
[10789]211               zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
[6140]212            END DO
213         END DO
214         DO jj = 2, jpjm1
215            DO ji = fs_2, fs_jpim1   ! vector opt.
[10789]216               zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( uu(ji,jj,jk,ktlev2) + uu(ji,jj,jk-1,ktlev2) )
217               zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vv(ji,jj,jk,ktlev2) + vv(ji,jj,jk-1,ktlev2) )
[6140]218            END DO
219         END DO
[643]220      END DO
[6140]221      DO jk = 1, jpkm1                          ! divergence of vertical momentum flux divergence
222         DO jj = 2, jpjm1
[643]223            DO ji = fs_2, fs_jpim1   ! vector opt.
[10789]224               pu_rhs(ji,jj,jk) =  pu_rhs(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev2)
225               pv_rhs(ji,jj,jk) =  pv_rhs(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev2)
[643]226            END DO
227         END DO
228      END DO
[1566]229      !
[6140]230      IF( l_trddyn ) THEN                       ! save the vertical advection trend for diagnostic
[10789]231         zfu_t(:,:,:) = pu_rhs(:,:,:) - zfu_t(:,:,:)
232         zfv_t(:,:,:) = pv_rhs(:,:,:) - zfv_t(:,:,:)
[5062]233         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
[1129]234      ENDIF
[6140]235      !                                         ! Control print
[10802]236      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
[1129]237         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]238      !
[643]239   END SUBROUTINE dyn_adv_ubs
240
241   !!==============================================================================
242END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.