Ticket #361: traadv_tvd.F90

File traadv_tvd.F90, 22.5 KB (added by charris, 11 years ago)

Fixed version of tvd advection code for poleward transport diagnostics

Line 
1MODULE traadv_tvd
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_tvd  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :       !  95-12  (L. Mortier)  Original code
7   !!                 !  00-01  (H. Loukos)  adapted to ORCA
8   !!                 !  00-10  (MA Foujols E.Kestenare)  include file not routine
9   !!                 !  00-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes
10   !!                 !  01-07  (E. Durand G. Madec)  adaptation to ORCA config
11   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
12   !!            9.0  !  04-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
13   !!            9.0  !  08-04  (S. Cravatte) add the i-, j- & k- trends computation
14   !!            " "  !  05-11  (V. Garnier) Surface pressure gradient organization
15   !!----------------------------------------------------------------------
16
17
18   !!----------------------------------------------------------------------
19   !!   tra_adv_tvd  : update the tracer trend with the horizontal
20   !!                  and vertical advection trends using a TVD scheme
21   !!   nonosc       : compute monotonic tracer fluxes by a nonoscillatory
22   !!                  algorithm
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean dynamics and active tracers
25   USE dom_oce         ! ocean space and time domain
26   USE trdmod          ! ocean active tracers trends
27   USE trdmod_oce      ! ocean variables trends
28   USE in_out_manager  ! I/O manager
29   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
30   USE trabbl          ! Advective term of BBL
31   USE lib_mpp
32   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
33   USE diaptr          ! poleward transport diagnostics
34   USE prtctl          ! Print control
35
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   tra_adv_tvd    ! routine called by step.F90
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !!   OPA 9.0 , LOCEAN-IPSL (2006)
47   !! $Id$
48   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50
51CONTAINS
52
53   SUBROUTINE tra_adv_tvd( kt, pun, pvn, pwn )
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_adv_tvd  ***
56      !!
57      !! **  Purpose :   Compute the now trend due to total advection of
58      !!       tracers and add it to the general trend of tracer equations
59      !!
60      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
61      !!       corrected flux (monotonic correction)
62      !!       note: - this advection scheme needs a leap-frog time scheme
63      !!
64      !! ** Action : - update (ta,sa) with the now advective tracer trends
65      !!             - save the trends in (ztrdt,ztrds) ('key_trdtra')
66      !!----------------------------------------------------------------------
67      USE oce              , ztrdt => ua   ! use ua as workspace
68      USE oce              , ztrds => va   ! use va as workspace
69      !!
70      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
71      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component
72      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component
73      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component
74      !!
75      INTEGER  ::   ji, jj, jk              ! dummy loop indices
76      REAL(wp) ::                        &  ! temporary scalar
77         ztat, zsat,                     &  !    "         "   
78         z_hdivn_x, z_hdivn_y, z_hdivn
79      REAL(wp) ::   &
80         z2dtt, zbtr, zeu, zev,          &  ! temporary scalar
81         zew, z2, zbtr1,                 &  ! temporary scalar
82         zfp_ui, zfp_vj, zfp_wk,         &  !    "         "
83         zfm_ui, zfm_vj, zfm_wk             !    "         "
84      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztu, ztv, ztw   ! temporary workspace
85      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zsi, zsu, zsv, zsw   !    "           "
86      !!----------------------------------------------------------------------
87
88      zti(:,:,:) = 0.e0   ;   zsi(:,:,:) = 0.e0
89
90      IF( kt == nit000 .AND. lwp ) THEN
91         WRITE(numout,*)
92         WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme'
93         WRITE(numout,*) '~~~~~~~~~~~'
94      ENDIF
95
96      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
97      ELSE                                        ;    z2 = 2.
98      ENDIF
99
100      ! 1. Bottom value : flux set to zero
101      ! ---------------
102      ztu(:,:,jpk) = 0.e0   ;   zsu(:,:,jpk) = 0.e0
103      ztv(:,:,jpk) = 0.e0   ;   zsv(:,:,jpk) = 0.e0
104      ztw(:,:,jpk) = 0.e0   ;   zsw(:,:,jpk) = 0.e0
105      zti(:,:,jpk) = 0.e0   ;   zsi(:,:,jpk) = 0.e0
106
107
108      ! 2. upstream advection with initial mass fluxes & intermediate update
109      ! --------------------------------------------------------------------
110      ! upstream tracer flux in the i and j direction
111      DO jk = 1, jpkm1
112         DO jj = 1, jpjm1
113            DO ji = 1, fs_jpim1   ! vector opt.
114               zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
115               zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
116               ! upstream scheme
117               zfp_ui = zeu + ABS( zeu )
118               zfm_ui = zeu - ABS( zeu )
119               zfp_vj = zev + ABS( zev )
120               zfm_vj = zev - ABS( zev )
121               ztu(ji,jj,jk) = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj  ,jk)
122               ztv(ji,jj,jk) = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji  ,jj+1,jk)
123               zsu(ji,jj,jk) = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk)
124               zsv(ji,jj,jk) = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk)
125            END DO
126         END DO
127      END DO
128
129      ! upstream tracer flux in the k direction
130      ! Surface value
131      IF( lk_vvl ) THEN
132         ! variable volume : flux set to zero
133         ztw(:,:,1) = 0.e0
134         zsw(:,:,1) = 0.e0
135      ELSE
136         ! free surface-constant volume
137         DO jj = 1, jpj
138            DO ji = 1, jpi
139               zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1)
140               ztw(ji,jj,1) = zew * tb(ji,jj,1)
141               zsw(ji,jj,1) = zew * sb(ji,jj,1)
142            END DO
143         END DO
144      ENDIF
145
146      ! Interior value
147      DO jk = 2, jpkm1
148         DO jj = 1, jpj
149            DO ji = 1, jpi
150               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
151               zfp_wk = zew + ABS( zew )
152               zfm_wk = zew - ABS( zew )
153               ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1)
154               zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1)
155            END DO
156         END DO
157      END DO
158
159      ! total advective trend
160      DO jk = 1, jpkm1
161         z2dtt = z2 * rdttra(jk)
162         DO jj = 2, jpjm1
163            DO ji = fs_2, fs_jpim1   ! vector opt.
164               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
165               ! total intermediate advective trends
166               ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )   &
167                  &     + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )   &
168                  &     + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr
169               zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  )   & 
170                  &     + zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  )   &
171                  &     + zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1) ) * zbtr
172               ! update and guess with monotonic sheme
173               ta(ji,jj,jk) =  ta(ji,jj,jk) + ztat
174               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsat
175               zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ztat ) * tmask(ji,jj,jk)
176               zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * zsat ) * tmask(ji,jj,jk)
177            END DO
178         END DO
179      END DO
180
181      ! "zonal" mean advective heat and salt transport
182      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
183         pht_adv(:) = ptr_vj( ztv(:,:,:) )
184         pst_adv(:) = ptr_vj( zsv(:,:,:) )
185      ENDIF
186
187      ! Save the intermediate i / j / k advective trends for diagnostics
188      ! -------------------------------------------------------------------
189      ! Warning : We should use zun instead of un in the computations below, but we
190      ! also use hdivn which is computed with un, vn (check ???). So we use un, vn
191      ! for consistency. Results are therefore approximate with key_trabbl_adv.
192
193      IF( l_trdtra ) THEN
194         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
195         !
196         ! T/S ZONAL advection trends
197         DO jk = 1, jpkm1
198            DO jj = 2, jpjm1
199               DO ji = fs_2, fs_jpim1   ! vector opt.
200                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
201                  ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr
202                  ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr
203               END DO
204            END DO
205         END DO
206         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)    ! save the trends
207         !
208         ! T/S MERIDIONAL advection trends
209         DO jk = 1, jpkm1
210            DO jj = 2, jpjm1
211               DO ji = fs_2, fs_jpim1   ! vector opt.
212                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
213                  ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr
214                  ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr
215               END DO
216            END DO
217         END DO
218         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)     ! save the trends
219         !
220         ! T/S VERTICAL advection trends
221         DO jk = 1, jpkm1
222            DO jj = 2, jpjm1
223               DO ji = fs_2, fs_jpim1   ! vector opt.         
224                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
225                  ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
226                  ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr
227               END DO
228            END DO
229         END DO
230         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)     ! save the trends
231         !
232      ENDIF
233
234      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
235      CALL lbc_lnk( zti, 'T', 1. )
236      CALL lbc_lnk( zsi, 'T', 1. )
237
238
239      ! 3. antidiffusive flux : high order minus low order
240      ! --------------------------------------------------
241      ! antidiffusive flux on i and j
242      DO jk = 1, jpkm1
243         DO jj = 1, jpjm1
244            DO ji = 1, fs_jpim1   ! vector opt.
245               zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
246               zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
247               ztu(ji,jj,jk) = zeu * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) ) - ztu(ji,jj,jk)
248               zsu(ji,jj,jk) = zeu * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) ) - zsu(ji,jj,jk)
249               ztv(ji,jj,jk) = zev * ( tn(ji,jj,jk) + tn(ji,jj+1,jk) ) - ztv(ji,jj,jk)
250               zsv(ji,jj,jk) = zev * ( sn(ji,jj,jk) + sn(ji,jj+1,jk) ) - zsv(ji,jj,jk)
251            END DO
252         END DO
253      END DO
254     
255      ! antidiffusive flux on k
256      ! Surface value
257      ztw(:,:,1) = 0.e0
258      zsw(:,:,1) = 0.e0
259
260      ! Interior value
261      DO jk = 2, jpkm1
262         DO jj = 1, jpj
263            DO ji = 1, jpi
264               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
265               ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
266               zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk)
267            END DO
268         END DO
269      END DO
270
271      ! Lateral bondary conditions
272      CALL lbc_lnk( ztu, 'U', -1. )   ;   CALL lbc_lnk( zsu, 'U', -1. )
273      CALL lbc_lnk( ztv, 'V', -1. )   ;   CALL lbc_lnk( zsv, 'V', -1. )
274      CALL lbc_lnk( ztw, 'W',  1. )   ;   CALL lbc_lnk( zsw, 'W',  1. )
275
276      ! 4. monotonicity algorithm
277      ! -------------------------
278      CALL nonosc( tb, ztu, ztv, ztw, zti, z2 )
279      CALL nonosc( sb, zsu, zsv, zsw, zsi, z2 )
280
281
282      ! 5. final trend with corrected fluxes
283      ! ------------------------------------
284      DO jk = 1, jpkm1
285         DO jj = 2, jpjm1
286            DO ji = fs_2, fs_jpim1   ! vector opt. 
287               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
288               ! total advective trends
289               ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )   &
290                  &     + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )   &
291                  &     + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr
292               zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  )   &
293                  &     + zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  )   &
294                  &     + zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1) ) * zbtr
295               ! add them to the general tracer trends
296               ta(ji,jj,jk) = ta(ji,jj,jk) + ztat
297               sa(ji,jj,jk) = sa(ji,jj,jk) + zsat
298            END DO
299         END DO
300      END DO
301
302
303      ! Save the advective trends for diagnostics
304      ! --------------------------------------------
305
306      IF( l_trdtra ) THEN
307         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
308         !
309         ! T/S ZONAL advection trends
310         DO jk = 1, jpkm1
311            DO jj = 2, jpjm1
312               DO ji = fs_2, fs_jpim1   ! vector opt.
313                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
314                  !   N.B. This computation is not valid along OBCs (if any)
315                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
316                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
317                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
318                  !-- Compute T/S zonal advection trends
319                  ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x
320                  ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x
321               END DO
322            END DO
323         END DO
324         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt, cnbpas='bis')   ! <<< ADD TO PREVIOUSLY COMPUTED
325         !
326         ! T/S MERIDIONAL advection trends
327         DO jk = 1, jpkm1
328            DO jj = 2, jpjm1
329               DO ji = fs_2, fs_jpim1   ! vector opt.
330                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
331                  !   N.B. This computation is not valid along OBCs (if any)
332                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
333                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
334                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
335                  !-- Compute T/S meridional advection trends
336                  ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y         
337                  ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y         
338               END DO
339            END DO
340         END DO
341         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt, cnbpas='bis')   ! <<< ADD TO PREVIOUSLY COMPUTED
342         !
343         ! T/S VERTICAL advection trends
344         DO jk = 1, jpkm1
345            DO jj = 2, jpjm1
346               DO ji = fs_2, fs_jpim1   ! vector opt.
347                  zbtr1     = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
348#if defined key_zco
349                  zbtr      = zbtr1
350                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
351                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
352#else
353                  zbtr      = zbtr1 / fse3t(ji,jj,jk)
354                  z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk)
355                  z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk)
356#endif
357                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
358                  zbtr      = zbtr1 / fse3t(ji,jj,jk)
359                  ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr - tn(ji,jj,jk) * z_hdivn
360                  ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr - sn(ji,jj,jk) * z_hdivn
361               END DO
362            END DO
363         END DO
364         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt, cnbpas='bis')   ! <<< ADD TO PREVIOUSLY COMPUTED
365         !
366      ENDIF
367
368      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' tvd adv  - Ta: ', mask1=tmask,   &
369         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
370
371      ! "zonal" mean advective heat and salt transport
372      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
373         pht_adv(:) = ptr_vj( ztv(:,:,:) ) + pht_adv(:)
374         pst_adv(:) = ptr_vj( zsv(:,:,:) ) + pst_adv(:)
375      ENDIF
376      !
377   END SUBROUTINE tra_adv_tvd
378
379
380   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, prdt )
381      !!---------------------------------------------------------------------
382      !!                    ***  ROUTINE nonosc  ***
383      !!     
384      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
385      !!       scheme and the before field by a nonoscillatory algorithm
386      !!
387      !! **  Method  :   ... ???
388      !!       warning : pbef and paft must be masked, but the boundaries
389      !!       conditions on the fluxes are not necessary zalezak (1979)
390      !!       drange (1995) multi-dimensional forward-in-time and upstream-
391      !!       in-space based differencing for fluid
392      !!----------------------------------------------------------------------
393      REAL(wp), INTENT( in ) ::   prdt   ! ???
394      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( inout ) ::   &
395         pbef,                            & ! before field
396         paft,                            & ! after field
397         paa,                             & ! monotonic flux in the i direction
398         pbb,                             & ! monotonic flux in the j direction
399         pcc                                ! monotonic flux in the k direction
400      !!
401      INTEGER ::   ji, jj, jk               ! dummy loop indices
402      INTEGER ::   ikm1
403      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
404      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbup, zbdo
405      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
406      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv
407      REAL(wp) ::   zup, zdo
408      !!----------------------------------------------------------------------
409
410      zbig = 1.e+40
411      zrtrn = 1.e-15
412      zbetup(:,:,jpk) = 0.e0   ;   zbetdo(:,:,jpk) = 0.e0
413
414
415      ! Search local extrema
416      ! --------------------
417      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
418      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   &
419         &        paft * tmask - zbig * ( 1.e0 - tmask )  )
420      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   &
421         &        paft * tmask + zbig * ( 1.e0 - tmask )  )
422
423      DO jk = 1, jpkm1
424         ikm1 = MAX(jk-1,1)
425         z2dtt = prdt * rdttra(jk)
426         DO jj = 2, jpjm1
427            DO ji = fs_2, fs_jpim1   ! vector opt.
428
429               ! search maximum in neighbourhood
430               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
431                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
432                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   &
433                  &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  )
434
435               ! search minimum in neighbourhood
436               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
437                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
438                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   &
439                  &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  )
440
441               ! positive part of the flux
442               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
443                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
444                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
445
446               ! negative part of the flux
447               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
448                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
449                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
450
451               ! up & down beta terms
452               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
453               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
454               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
455            END DO
456         END DO
457      END DO
458
459      ! lateral boundary condition on zbetup & zbetdo   (unchanged sign)
460      CALL lbc_lnk( zbetup, 'T', 1. )
461      CALL lbc_lnk( zbetdo, 'T', 1. )
462
463
464      ! 3. monotonic flux in the i & j direction (paa & pbb)
465      ! ----------------------------------------
466      DO jk = 1, jpkm1
467         DO jj = 2, jpjm1
468            DO ji = fs_2, fs_jpim1   ! vector opt.
469               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
470               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
471               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
472               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
473
474               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
475               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
476               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
477               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
478
479      ! monotonic flux in the k direction, i.e. pcc
480      ! -------------------------------------------
481               za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
482               zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
483               zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
484               pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb )
485            END DO
486         END DO
487      END DO
488
489      ! lateral boundary condition on paa, pbb, pcc
490      CALL lbc_lnk( paa, 'U', -1. )      ! changed sign
491      CALL lbc_lnk( pbb, 'V', -1. )      ! changed sign
492      !
493   END SUBROUTINE nonosc
494
495   !!======================================================================
496END MODULE traadv_tvd