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

source: branches/dev_001_GM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 @ 2291

Last change on this file since 2291 was 790, checked in by gm, 16 years ago

dev_001_GM - complete theprevious comit with omitted routines

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 12.3 KB
Line 
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   !!======================================================================
7   !! History :  9.0  !  06-08  (R. Benshila, L. Debreu)  Original code
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   dyn_adv_ubs   : flux form momentum advection using    (ln_dynadv=T)
12   !!                   an 3rd order Upstream Biased Scheme or Quick scheme
13   !!                   combined with 2nd or 4th order finite differences
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers
16   USE dom_oce        ! ocean space and time domain
17   USE dynspg_oce     ! surface pressure gradient
18   USE in_out_manager ! I/O manager
19   USE dynspg_rl      ! I/O manager
20   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
21
22   IMPLICIT NONE
23   PRIVATE
24
25   REAL(wp), PARAMETER :: gamma1 = 1._wp/4._wp  ! =1/4 quick      ; =1/3  3rd order UBS
26   REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp  ! =0   2nd order  ; =1/8  4th order centred
27
28   !! * Routine accessibility
29   PUBLIC dyn_adv_ubs                            ! routine called by step.F90
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !!   OPA 9.0 , LODYC-IPSL  (2006)
36   !! $Header$
37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE dyn_adv_ubs( kt )
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE dyn_adv_ubs  ***
45      !!
46      !! ** Purpose :   Compute the now momentum advection trend in flux form
47      !!      and the general trend of the momentum equation.
48      !!
49      !! ** Method  :   The scheme is the one implemeted in ROMS. It depends
50      !!      on two parameter gamma1 and gamma2. The former control the
51      !!      upstream baised part of the scheme and the later the centred
52      !!      part:     gamma1 = 0    pure centered  (no diffusive part)
53      !!                       = 1/4  Quick scheme
54      !!                       = 1/3  3rd order Upstream biased scheme
55      !!                gamma2 = 0    2nd order finite differencing
56      !!                       = 1/8  4th order finite differencing
57      !!      For stability reasons, the first term of the fluxes which cor-
58      !!      responds to a second order centered scheme is evaluated using 
59      !!      the now velocity (centered in time) while the second term which 
60      !!      is the diffusive part of the scheme, is evaluated using the
61      !!      before velocity (forward in time).
62      !!      Default value (hard coded in the begining of the module) are
63      !!      gamma1=1/4 and gamma2=1/8.
64      !!
65      !! ** Action : - update (ua,va) with the 3D advective momentum trends
66      !!
67      !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
68      !!----------------------------------------------------------------------
69      USE oce, ONLY:   zfu => ta,   & ! use ta as 3D workspace
70                       zfv => sa      ! use sa as 3D workspace
71
72      INTEGER, INTENT(in) ::   kt     ! ocean time-step index
73
74      INTEGER  ::   ji, jj, jk                      ! dummy loop indices
75      REAL(wp) ::   zua, zva, zbu, zbv              ! temporary scalars
76      REAL(wp) ::   zui, zvj, zfuj, zfvi, zl_u, zl_v! temporary scalars
77      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfu_t, zfu_f     ! temporary workspace
78      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfv_t, zfv_f     !    "           "
79      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zfw, zfu_uw, zfv_vw
80      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlu_uu, zlu_uv   ! temporary workspace
81      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   zlv_vv, zlv_vu   ! temporary workspace
82      !!----------------------------------------------------------------------
83
84      IF( kt == nit000 ) THEN
85         IF(lwp) WRITE(numout,*)
86         IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
87         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
88      ENDIF
89      zfu_t(:,:,:) = 0.e0
90      zfv_t(:,:,:) = 0.e0
91      zfu_f(:,:,:) = 0.e0
92      zfv_f(:,:,:) = 0.e0
93     
94      zlu_uu(:,:,:,:) = 0.e0 
95      zlv_vv(:,:,:,:) = 0.e0 
96      zlu_uv(:,:,:,:) = 0.e0 
97      zlv_vu(:,:,:,:) = 0.e0 
98
99     
100      !                                                ! ===============
101      DO jk = 1, jpkm1                                 ! Horizontal slab
102         !                                             ! ===============
103         ! Laplacian
104         ! ---------
105         zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
106         zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
107           
108         DO jj = 2, jpjm1
109            DO ji = fs_2, fs_jpim1   ! vector opt.
110               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)
111               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) 
112               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)
113               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)
114               
115               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)
116               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)
117               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)
118               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)
119            END DO
120         END DO
121         !                                             ! ===============
122      END DO                                           !   End of slab
123      !                                                ! ===============
124
125      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.) 
126      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.) 
127      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.) 
128      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.) 
129
130      CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) 
131      CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) 
132      CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) 
133      CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 
134
135      ! I. Horizontal advection
136      ! -----------------------
137      !                                                ! ===============
138      DO jk = 1, jpkm1                                 ! Horizontal slab
139         !                                             ! ===============
140         ! horizontal volume fluxes
141         zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
142         zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
143
144         ! horizontal momentum fluxes at T- and F-point
145         DO jj = 1, jpjm1
146            DO ji = 1, fs_jpim1   ! vector opt.
147               zui = ( un(ji,jj,jk) + un(ji+1,jj  ,jk) )
148               zvj = ( vn(ji,jj,jk) + vn(ji  ,jj+1,jk) )
149
150               IF (zui > 0) THEN   ;   zl_u = zlu_uu(ji  ,jj,jk,1)
151               ELSE                ;   zl_u = zlu_uu(ji+1,jj,jk,1)
152               ENDIF
153               IF (zvj > 0) THEN   ;   zl_v = zlv_vv(ji,jj  ,jk,1)
154               ELSE                ;   zl_v = zlv_vv(ji,jj+1,jk,1)
155               ENDIF
156
157               zfu_t(ji+1,jj  ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj  ,jk)                               &
158                  &                    - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj  ,jk,2) )  )   &
159                  &                * ( zui - gamma1 * zl_u)
160               zfv_t(ji  ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji  ,jj+1,jk)                               &
161                  &                    - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji  ,jj+1,jk,2) )  )   &
162                  &                * ( zvj - gamma1 * zl_v)
163
164               zfuj = ( zfu(ji,jj,jk) + zfu(ji  ,jj+1,jk) )
165               zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj  ,jk) )
166               IF (zfuj > 0) THEN   ;    zl_v = zlv_vu( ji  ,jj  ,jk,1)
167               ELSE                 ;    zl_v = zlv_vu( ji+1,jj,jk,1)
168               ENDIF
169               IF (zfvi > 0) THEN   ;    zl_u = zlu_uv( ji,jj  ,jk,1)
170               ELSE                 ;    zl_u = zlu_uv( ji,jj+1,jk,1)
171               ENDIF
172
173               zfv_f(ji  ,jj  ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj  ,jk,2) )  )   &
174                  &                * ( un(ji,jj,jk) + un(ji  ,jj+1,jk) - gamma1 * zl_u )
175               zfu_f(ji  ,jj  ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji  ,jj+1,jk,2) )  )   &
176                  &                * ( vn(ji,jj,jk) + vn(ji+1,jj  ,jk) - gamma1 * zl_v )
177            END DO
178         END DO
179
180         ! divergence of horizontal momentum fluxes
181         DO jj = 2, jpjm1
182            DO ji = fs_2, fs_jpim1   ! vector opt.
183               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
184               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
185               ! horizontal advective trends
186               zua = - (  zfu_t(ji+1,jj  ,jk) - zfu_t(ji  ,jj  ,jk)   &
187                  &     + zfv_f(ji  ,jj  ,jk) - zfv_f(ji  ,jj-1,jk)  ) / zbu
188               zva = - (  zfu_f(ji  ,jj  ,jk) - zfu_f(ji-1,jj  ,jk)   &
189                  &     + zfv_t(ji  ,jj+1,jk) - zfv_t(ji  ,jj  ,jk)  ) / zbv
190               ! add it to the general tracer trends
191               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
192               va(ji,jj,jk) = va(ji,jj,jk) + zva
193            END DO
194         END DO
195         !                                             ! ===============
196      END DO                                           !   End of slab
197      !                                                ! ===============
198
199!!gm  momentum trend diagnostics is missing
200
201      ! II. Vertical advection
202      ! ----------------------
203
204      ! Second order centered tracer flux at w-point
205      DO jk = 1, jpkm1
206         ! Vertical volume fluxes                   
207         zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk)
208         ! Vertical advective fluxes                   
209         IF( jk == 1 ) THEN
210            zfu_uw(:,:,jpk) = 0.e0         ! Bottom value : flux set to zero
211            zfv_vw(:,:,jpk) = 0.e0
212            !                              ! Surface value
213            IF( lk_dynspg_rl ) THEN             ! rigid lid : flux set to zero
214               zfu_uw(:,:, 1 ) = 0.e0   
215               zfv_vw(:,:, 1 ) = 0.e0
216            ELSE                                ! free surface-constant volume
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
224         ELSE
225            !                              ! interior fluxes
226            DO jj = 2, jpjm1
227               DO ji = fs_2, fs_jpim1   ! vector opt.
228                  zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj  ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )
229                  zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji  ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )
230               END DO
231            END DO
232         ENDIF
233      END DO
234
235      ! momentum flux divergence at u-, v-points added to the general trend
236      DO jk = 1, jpkm1
237         DO jj = 2, jpjm1 
238            DO ji = fs_2, fs_jpim1   ! vector opt.
239               zua = - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) )    &
240                  &  / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
241               zva = - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) )    &
242                  &  / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
243               ! add it to the general tracer trends
244               ua(ji,jj,jk) =  ua(ji,jj,jk) + zua
245               va(ji,jj,jk) =  va(ji,jj,jk) + zva
246            END DO
247         END DO
248      END DO
249
250   END SUBROUTINE dyn_adv_ubs
251
252   !!==============================================================================
253END MODULE dynadv_ubs
Note: See TracBrowser for help on using the repository browser.