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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 2865

Last change on this file since 2865 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk
  • Property svn:keywords set to Id
File size: 14.1 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
[1129]18   USE trdmod         ! ocean dynamics trends
19   USE trdmod_oce     ! ocean variables trends
[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
[643]24
25   IMPLICIT NONE
26   PRIVATE
27
28   REAL(wp), PARAMETER :: gamma1 = 1._wp/4._wp  ! =1/4 quick      ; =1/3  3rd order UBS
29   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred
30
[1566]31   PUBLIC   dyn_adv_ubs   ! routine called by step.F90
[643]32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
[2715]37   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
[1152]38   !! $Id$
[2715]39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[643]40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE dyn_adv_ubs( kt )
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
57      !!                       = 1/8  4th order finite differencing
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
64      !!      gamma1=1/4 and gamma2=1/8.
65      !!
[1566]66      !! ** Action : - (ua,va) updated with the 3D advective momentum trends
[643]67      !!
68      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
69      !!----------------------------------------------------------------------
[2715]70      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
71      USE oce     , ONLY:   zfu    => ta       , zfv    => sa      ! (ta,sa) used as 3D workspace
72      USE wrk_nemo, ONLY:   zfu_t  => wrk_3d_1 , zfv_t  =>wrk_3d_4 , zfu_uw =>wrk_3d_6   ! 3D workspace
73      USE wrk_nemo, ONLY:   zfu_f  => wrk_3d_2 , zfv_f  =>wrk_3d_5 , zfv_vw =>wrk_3d_7
74      USE wrk_nemo, ONLY:   zfw    => wrk_3d_3
75      USE wrk_nemo, ONLY:   zlu_uu => wrk_4d_1 , zlv_vv=>wrk_4d_3   ! 4D workspace
76      USE wrk_nemo, ONLY:   zlu_uv => wrk_4d_2 , zlv_vu=>wrk_4d_4
77      !
[643]78      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
[2715]79      !
[1566]80      INTEGER  ::   ji, jj, jk            ! dummy loop indices
81      REAL(wp) ::   zbu, zbv    ! temporary scalars
82      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v   ! temporary scalars
[643]83      !!----------------------------------------------------------------------
84
85      IF( kt == nit000 ) THEN
86         IF(lwp) WRITE(numout,*)
87         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
88         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
89      ENDIF
[2715]90
91      ! Check that required workspace arrays are not already in use
92      IF( wrk_in_use(3, 1,2,3,4,5,6,7) .OR. wrk_in_use(4, 1,2,3,4) ) THEN
93         CALL ctl_stop('dyn_adv_ubs: requested workspace array unavailable')   ;   RETURN
94      ENDIF
95
96      zfu_t(:,:,:) = 0._wp
97      zfv_t(:,:,:) = 0._wp
98      zfu_f(:,:,:) = 0._wp
99      zfv_f(:,:,:) = 0._wp
[1566]100      !
[2715]101      zlu_uu(:,:,:,:) = 0._wp
102      zlv_vv(:,:,:,:) = 0._wp 
103      zlu_uv(:,:,:,:) = 0._wp 
104      zlv_vu(:,:,:,:) = 0._wp 
[643]105
[1129]106      IF( l_trddyn ) THEN           ! Save ua and va trends
107         zfu_uw(:,:,:) = ua(:,:,:)
108         zfv_vw(:,:,:) = va(:,:,:)
109      ENDIF
110
[1566]111      !                                      ! =========================== !
112      DO jk = 1, jpkm1                       !  Laplacian of the velocity  !
113         !                                   ! =========================== !
114         !                                         ! horizontal volume fluxes
[643]115         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
116         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]117         !           
118         DO jj = 2, jpjm1                          ! laplacian
[643]119            DO ji = fs_2, fs_jpim1   ! vector opt.
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)-2.*ub (ji,jj,jk)+ub (ji,jj-1,jk) ) * umask(ji,jj,jk)
123               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]124               !
[643]125               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)
126               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)
127               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)
128               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)
129            END DO
130         END DO
[1566]131      END DO
132!!gm BUG !!!  just below this should be +1 in all the communications
[2715]133!      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.)
134!      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.)
135!      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.)
136!      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.)
137!
[1566]138!!gm corrected:
139      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. )
140      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. )
141      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. )
142      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 
143!!gm end
144     
145      !                                      ! ====================== !
146      !                                      !  Horizontal advection  !
147      DO jk = 1, jpkm1                       ! ====================== !
148         !                                         ! horizontal volume fluxes
[643]149         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
150         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
[1566]151         !
152         DO jj = 1, jpjm1                          ! horizontal momentum fluxes at T- and F-point
[643]153            DO ji = 1, fs_jpim1   ! vector opt.
154               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
155               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
[1566]156               !
[643]157               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
158               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
159               ENDIF
160               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
161               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
162               ENDIF
[1566]163               !
[643]164               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
165                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
166                  &                * ( zui - gamma1 * zl_u)
167               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
168                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
169                  &                * ( zvj - gamma1 * zl_v)
[1566]170               !
[643]171               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
172               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
173               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
174               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
175               ENDIF
176               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
177               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
178               ENDIF
[1566]179               !
[643]180               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
181                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
182               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
183                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
184            END DO
185         END DO
[1566]186         DO jj = 2, jpjm1                          ! divergence of horizontal momentum fluxes
[643]187            DO ji = fs_2, fs_jpim1   ! vector opt.
188               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
189               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
[1566]190               !
191               ua(ji,jj,jk) = ua(ji,jj,jk) - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)    &
192                  &                           + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
193               va(ji,jj,jk) = va(ji,jj,jk) - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)    &
194                  &                           + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
[643]195            END DO
196         END DO
[1566]197      END DO
198      IF( l_trddyn ) THEN                          ! save the horizontal advection trend for diagnostic
[1129]199         zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
200         zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
201         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
202         zfu_t(:,:,:) = ua(:,:,:)
203         zfv_t(:,:,:) = va(:,:,:)
204      ENDIF
205
[1566]206      !                                      ! ==================== !
207      !                                      !  Vertical advection  !
208      DO jk = 1, jpkm1                       ! ==================== !
209         !                                         ! Vertical volume fluxesÊ
[643]210         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
[1566]211         !
212         IF( jk == 1 ) THEN                        ! surface/bottom advective fluxes                   
213            zfu_uw(:,:,jpk) = 0.e0                      ! Bottom  value : flux set to zero
[643]214            zfv_vw(:,:,jpk) = 0.e0
[1566]215            !                                           ! Surface value :
216            IF( lk_vvl ) THEN                                ! variable volume : flux set to zero
[643]217               zfu_uw(:,:, 1 ) = 0.e0   
218               zfv_vw(:,:, 1 ) = 0.e0
[1566]219            ELSE                                             ! constant volume : advection through the surface
[643]220               DO jj = 2, jpjm1
221                  DO ji = fs_2, fs_jpim1
222                     zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj  ,1) ) * un(ji,jj,1)
223                     zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji  ,jj+1,1) ) * vn(ji,jj,1)
224                  END DO
225               END DO
226            ENDIF
[1566]227         ELSE                                      ! interior fluxes
[643]228            DO jj = 2, jpjm1
229               DO ji = fs_2, fs_jpim1   ! vector opt.
230                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
231                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
232               END DO
233            END DO
234         ENDIF
235      END DO
[1566]236      DO jk = 1, jpkm1                             ! divergence of vertical momentum flux divergence
[643]237         DO jj = 2, jpjm1 
238            DO ji = fs_2, fs_jpim1   ! vector opt.
[1566]239               ua(ji,jj,jk) =  ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
[643]240                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
[1566]241               va(ji,jj,jk) =  va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
[643]242                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
243            END DO
244         END DO
245      END DO
[1566]246      !
247      IF( l_trddyn ) THEN                          ! save the vertical advection trend for diagnostic
[1129]248         zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
249         zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
250         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
251      ENDIF
[1566]252      !                                            ! Control print
[1129]253      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
254         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
[1566]255      !
[2715]256      IF( wrk_not_released(3, 1,2,3,4,5,6,7) .OR.   &
257          wrk_not_released(4, 1,2,3,4)       )   CALL ctl_stop('dyn_adv_ubs: failed to release workspace array')
258      !
[643]259   END SUBROUTINE dyn_adv_ubs
260
261   !!==============================================================================
262END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.