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/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 3318

Last change on this file since 3318 was 3318, checked in by gm, 12 years ago

Ediag branche: #927 split TRA/DYN trd computation
  • Property svn:keywords set to Id
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
[3318]18   USE trd_oce        ! trends: ocean variables
19   USE trddyn         ! trend manager: dynamics
[2715]20   USE in_out_manager ! I/O manager
[1129]21   USE prtctl         ! Print control
[2715]22   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
23   USE lib_mpp        ! MPP library
[3294]24   USE wrk_nemo        ! Memory Allocation
25   USE timing          ! Timing
[643]26
27   IMPLICIT NONE
28   PRIVATE
29
[3294]30   REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp  ! =1/4 quick      ; =1/3  3rd order UBS
[643]31   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred
32
[1566]33   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
[643]34
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
[2715]39   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[1152]40   !! $Id$
[2715]41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[643]42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE dyn_adv_ubs( kt )
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
59      !!                       = 1/8  4th order finite differencing
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
[3294]66      !!      gamma1=1/3 and gamma2=1/8.
[643]67      !!
[1566]68      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
[643]69      !!
70      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
71      !!----------------------------------------------------------------------
72      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
[2715]73      !
[1566]74      INTEGER  ::   ji, jj, jk            ! dummy loop indices
75      REAL(wp) ::   zbu, zbv    ! temporary scalars
76      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! temporary scalars
[3294]77      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu, zfv
78      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::  zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw
79      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zlu_uu, zlv_vv, zlu_uv, zlv_vu
[643]80      !!----------------------------------------------------------------------
[3294]81      !
82      IF( nn_timing == 1 )  CALL timing_start('dyn_adv_ubs')
83      !
84      CALL wrk_alloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
85      CALL wrk_alloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                               )
86      !
[643]87      IF( kt == nit000 ) THEN
88         IF(lwp) WRITE(numout,*)
89         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
90         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
91      ENDIF
[3294]92      !
[2715]93      zfu_t(:,:,:) = 0._wp
94      zfv_t(:,:,:) = 0._wp
95      zfu_f(:,:,:) = 0._wp
96      zfv_f(:,:,:) = 0._wp
[1566]97      !
[2715]98      zlu_uu(:,:,:,:) = 0._wp
99      zlv_vv(:,:,:,:) = 0._wp 
100      zlu_uv(:,:,:,:) = 0._wp 
101      zlv_vu(:,:,:,:) = 0._wp 
[643]102
[1129]103      IF( l_trddyn ) THEN           ! Save ua and va trends
104         zfu_uw(:,:,:) = ua(:,:,:)
105         zfv_vw(:,:,:) = va(:,:,:)
106      ENDIF
107
[1566]108      !                                      ! =========================== !
109      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
110         !                                   ! =========================== !
111         !                                         ! horizontal volume fluxes
[643]112         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
113         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]114         !           
115         DO jj = 2, jpjm1                          ! laplacian
[643]116            DO ji = fs_2, fs_jpim1   ! vector opt.
117               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)
118               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) 
119               zlu_uv(ji,jj,jk,1) = ( ub (ji,jj+1,jk)-2.*ub (ji,jj,jk)+ub (ji,jj-1,jk) ) * umask(ji,jj,jk)
120               zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj,jk)-2.*vb (ji,jj,jk)+vb (ji-1,jj,jk) ) * vmask(ji,jj,jk)
[2715]121               !
[643]122               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)
123               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)
124               zlu_uv(ji,jj,jk,2) = ( zfu(ji,jj+1,jk)-2.*zfu(ji,jj,jk)+zfu(ji,jj-1,jk) ) * umask(ji,jj,jk)
125               zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj,jk)-2.*zfv(ji,jj,jk)+zfv(ji-1,jj,jk) ) * vmask(ji,jj,jk)
126            END DO
127         END DO
[1566]128      END DO
129!!gm BUG !!!  just below this should be +1 in all the communications
[2715]130!      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.)
131!      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.)
132!      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.)
133!      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.)
134!
[1566]135!!gm corrected:
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!!gm end
141     
142      !                                      ! ====================== !
143      !                                      !  Horizontal advection  !
144      DO jk = 1, jpkm1                       ! ====================== !
145         !                                         ! horizontal volume fluxes
[643]146         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
147         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]148         !
149         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
[643]150            DO ji = 1, fs_jpim1   ! vector opt.
151               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
152               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
[1566]153               !
[643]154               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
155               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
156               ENDIF
157               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
158               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
159               ENDIF
[1566]160               !
[643]161               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
162                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
163                  &                * ( zui - gamma1 * zl_u)
164               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
165                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
166                  &                * ( zvj - gamma1 * zl_v)
[1566]167               !
[643]168               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
169               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
170               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
171               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
172               ENDIF
173               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
174               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
175               ENDIF
[1566]176               !
[643]177               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
178                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
179               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
180                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
181            END DO
182         END DO
[1566]183         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
[643]184            DO ji = fs_2, fs_jpim1   ! vector opt.
185               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
186               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
[1566]187               !
188               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
189                  &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
190               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    &
191                  &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
[643]192            END DO
193         END DO
[1566]194      END DO
195      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
[1129]196         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
197         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
[3318]198         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
[1129]199         zfu_t(:,:,:) = ua(:,:,:)
200         zfv_t(:,:,:) = va(:,:,:)
201      ENDIF
202
[1566]203      !                                      ! ==================== !
204      !                                      !  Vertical advection  !
205      DO jk = 1, jpkm1                       ! ==================== !
206         !                                         ! Vertical volume fluxesÊ
[643]207         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
[1566]208         !
209         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
210            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
[643]211            zfv_vw(:,:,jpk) = 0.e0
[1566]212            !                                           ! Surface value :
213            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
[643]214               zfu_uw(:,:, 1 ) = 0.e0   
215               zfv_vw(:,:, 1 ) = 0.e0
[1566]216            ELSE                                             ! constant volume : advection through the surface
[643]217               DO jj = 2, jpjm1
218                  DO ji = fs_2, fs_jpim1
219                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
220                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
221                  END DO
222               END DO
223            ENDIF
[1566]224         ELSE                                      ! interior fluxes
[643]225            DO jj = 2, jpjm1
226               DO ji = fs_2, fs_jpim1   ! vector opt.
227                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
228                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
229               END DO
230            END DO
231         ENDIF
232      END DO
[1566]233      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
[643]234         DO jj = 2, jpjm1 
235            DO ji = fs_2, fs_jpim1   ! vector opt.
[1566]236               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
[643]237                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
[1566]238               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
[643]239                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
240            END DO
241         END DO
242      END DO
[1566]243      !
244      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
[1129]245         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
246         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
[3318]247         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
[1129]248      ENDIF
[1566]249      !                                            ! Control print
[1129]250      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
251         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]252      !
[3294]253      CALL wrk_dealloc( jpi, jpj, jpk,       zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )
254      CALL wrk_dealloc( jpi, jpj, jpk, jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu                               )
[2715]255      !
[3294]256      IF( nn_timing == 1 )  CALL timing_stop('dyn_adv_ubs')
257      !
[643]258   END SUBROUTINE dyn_adv_ubs
259
260   !!==============================================================================
261END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.