source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 3432

Last change on this file since 3432 was 3432, checked in by trackstand2, 9 years ago

Merge branch 'ksection_partition'

  • Property svn:keywords set to Id
File size: 21.9 KB
Line 
1MODULE traadv_tvd
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_tvd  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  OPA  !  1995-12  (L. Mortier)  Original code
7   !!                 !  2000-01  (H. Loukos)  adapted to ORCA
8   !!                 !  2000-10  (MA Foujols E.Kestenare)  include file not routine
9   !!                 !  2000-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes
10   !!                 !  2001-07  (E. Durand G. Madec)  adaptation to ORCA config
11   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
12   !!    NEMO    1.0  !  2004-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
13   !!            2.0  !  2008-04  (S. Cravatte) add the i-, j- & k- trends computation
14   !!             -   !  2009-11  (V. Garnier) Surface pressure gradient organization
15   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
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_oce      ! tracers trends
27   USE trdtra          ! tracers trends
28   USE in_out_manager  ! I/O manager
29   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
30   USE lib_mpp         ! MPP library
31   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
32   USE diaptr          ! poleward transport diagnostics
33   USE trc_oce         ! share passive tracers/Ocean variables
34
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   tra_adv_tvd    ! routine called by step.F90
40
41   LOGICAL ::   l_trd   ! flag to compute trends
42
43   !! * Control permutation of array indices
44#  include "oce_ftrans.h90"
45#  include "dom_oce_ftrans.h90"
46#  include "trc_oce_ftrans.h90"
47
48   !! * Substitutions
49#  include "domzgr_substitute.h90"
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
53   !! $Id$
54   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE tra_adv_tvd ( kt, cdtype, p2dt, pun, pvn, pwn,      &
59      &                                       ptb, ptn, pta, kjpt )
60      !!----------------------------------------------------------------------
61      !!                  ***  ROUTINE tra_adv_tvd  ***
62      !!
63      !! **  Purpose :   Compute the now trend due to total advection of
64      !!       tracers and add it to the general trend of tracer equations
65      !!
66      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
67      !!       corrected flux (monotonic correction)
68      !!       note: - this advection scheme needs a leap-frog time scheme
69      !!
70      !! ** Action : - update (pta) with the now advective tracer trends
71      !!             - save the trends
72      !!----------------------------------------------------------------------
73      USE timing,   ONLY: timing_start, timing_stop
74      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
75      USE oce     , ONLY:   zwx => ua        , zwy => va          ! (ua,va) used as workspace
76      USE wrk_nemo, ONLY:   zwi => wrk_3d_12 , zwz => wrk_3d_13   ! 3D workspace
77
78      !! DCSE_NEMO: need additional directives for renamed module variables
79!FTRANS zwx zwy zwi zwz :I :I :z
80
81      !
82      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
83      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
84      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
85      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
86
87      !! DCSE_NEMO: This style defeats ftrans
88!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
89!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
90!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
91
92!FTRANS pun pvn pwn :I :I :z
93!FTRANS ptb ptn pta :I :I :z :
94      REAL(wp), INTENT(in   ) ::   pun(jpi,jpj,jpk)        ! ocean velocity component (u)
95      REAL(wp), INTENT(in   ) ::   pvn(jpi,jpj,jpk)        ! ocean velocity component (v)
96      REAL(wp), INTENT(in   ) ::   pwn(jpi,jpj,jpk)        ! ocean velocity component (w)
97      REAL(wp), INTENT(in   ) ::   ptb(jpi,jpj,jpk,kjpt)   ! tracer fields (before)
98      REAL(wp), INTENT(in   ) ::   ptn(jpi,jpj,jpk,kjpt)   ! tracer fields (now)
99      REAL(wp), INTENT(inout) ::   pta(jpi,jpj,jpk,kjpt)   ! tracer trend
100
101      !
102      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
103      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar
104      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
105      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
106      REAL(wp), DIMENSION (:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz
107!FTRANS ztrdx ztrdy ztrdz :I :I :z
108
109      !!----------------------------------------------------------------------
110
111      CALL timing_start('tra_adv_tvd')
112
113      IF( wrk_in_use(3, 12,13) ) THEN
114         CALL ctl_stop('tra_adv_tvd: requested workspace arrays unavailable')   ;   RETURN
115      ENDIF
116
117      IF( kt == nit000 )  THEN
118         IF(lwp) WRITE(numout,*)
119         IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
120         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
121         !
122         l_trd = .FALSE.
123         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
124      ENDIF
125      !
126      IF( l_trd )  THEN
127        ALLOCATE( ztrdx(jpi,jpj,jpk) )      ;      ztrdx(:,:,:) = 0.e0
128        ALLOCATE( ztrdy(jpi,jpj,jpk) )      ;      ztrdy(:,:,:) = 0.e0
129        ALLOCATE( ztrdz(jpi,jpj,jpk) )      ;      ztrdz(:,:,:) = 0.e0
130      END IF
131      !
132      zwi(:,:,:) = 0.e0
133      !
134      !                                                          ! ===========
135      DO jn = 1, kjpt                                            ! tracer loop
136         !                                                       ! ===========
137         ! 1. Bottom value : flux set to zero
138         ! ----------------------------------
139         zwx(:,:,jpk) = 0.e0    ;    zwz(:,:,jpk) = 0.e0
140         zwy(:,:,jpk) = 0.e0    ;    zwi(:,:,jpk) = 0.e0
141
142         ! 2. upstream advection with initial mass fluxes & intermediate update
143         ! --------------------------------------------------------------------
144         ! upstream tracer flux in the i and j direction
145         CALL timing_start('tvd_upstream')
146#if defined key_z_first
147         DO jj = 1, jpjm1
148            DO ji = 1, jpim1
149               DO jk = 1, jpkm1
150#else
151         DO jk = 1, jpkm1
152            DO jj = 1, jpjm1
153               DO ji = 1, fs_jpim1   ! vector opt.
154#endif
155                  ! upstream scheme
156                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
157                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
158                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
159                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
160                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
161                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
162               END DO
163            END DO
164         END DO
165         CALL timing_stop('tvd_upstream','section')
166
167         ! upstream tracer flux in the k direction
168         ! Surface value
169         CALL timing_start('tvd_upstreamk')
170         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable
171         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface
172         ENDIF
173         ! Interior value
174#if defined key_z_first
175         DO jj = 1, jpj
176            DO ji = 1, jpi
177               DO jk = 2, jpkm1
178#else
179         DO jk = 2, jpkm1
180            DO jj = 1, jpj
181               DO ji = 1, jpi
182#endif
183                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
184                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
185                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
186               END DO
187            END DO
188         END DO
189         CALL timing_stop('tvd_upstreamk','section')
190
191         ! total advective trend
192         CALL timing_start('tvd_tot')
193#if defined key_z_first
194         DO jj = 2, jpjm1
195            DO ji = 2, jpim1
196               DO jk = 1, jpkm1
197                  z2dtt = p2dt(jk)
198#else
199         DO jk = 1, jpkm1
200            z2dtt = p2dt(jk)
201            DO jj = 2, jpjm1
202               DO ji = fs_2, fs_jpim1   ! vector opt.
203#endif
204                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
205                  ! total intermediate advective trends
206                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
207                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
208                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
209                  ! update and guess with monotonic sheme
210                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra
211                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
212               END DO
213            END DO
214         END DO
215         CALL timing_stop('tvd_tot','section')
216
217         !                             ! Lateral boundary conditions on zwi  (unchanged sign)
218         CALL timing_start('tvd_lbc')
219         CALL lbc_lnk( zwi, 'T', 1. ) 
220         CALL timing_stop('tvd_lbc','section')
221
222         !                                 ! trend diagnostics (contribution of upstream fluxes)
223         IF( l_trd )  THEN 
224            ! store intermediate advective trends
225            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
226         END IF
227         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
228         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
229           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
230           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
231         ENDIF
232
233         ! 3. antidiffusive flux : high order minus low order
234         ! --------------------------------------------------
235         ! antidiffusive flux on i and j
236         CALL timing_start('tvd_antidiff')
237#if defined key_z_first
238         DO jj = 1, jpjm1
239            DO ji = 1, jpim1
240               DO jk = 1, jpkm1
241#else
242         DO jk = 1, jpkm1
243            DO jj = 1, jpjm1
244               DO ji = 1, fs_jpim1   ! vector opt.
245#endif
246                  zwx(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
247                  zwy(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
248               END DO
249            END DO
250         END DO
251         CALL timing_stop('tvd_antidiff','section')
252     
253         ! antidiffusive flux on k
254         CALL timing_start('tvd_antidiffk')
255#if defined key_z_first
256         DO jj = 1, jpj
257            DO ji = 1, jpi
258               zwz(ji,jj,1) = 0.e0   ! Surface value
259               DO jk = 2, jpkm1
260#else
261         zwz(:,:,1) = 0.e0           ! Surface value
262         !
263         DO jk = 2, jpkm1            ! Interior value
264            DO jj = 1, jpj
265               DO ji = 1, jpi
266#endif
267                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk)
268               END DO
269            END DO
270         END DO
271         CALL timing_stop('tvd_antidiffk','section')
272
273         CALL timing_start('tvd_lbc')
274         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
275         CALL lbc_lnk( zwz, 'W',  1. )
276         CALL timing_stop('tvd_lbc','section')
277
278         ! 4. monotonicity algorithm
279         ! -------------------------
280         CALL timing_start('tvd_nonosc')
281         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
282         CALL timing_stop('tvd_nonosc','section')
283
284
285         ! 5. final trend with corrected fluxes
286         ! ------------------------------------
287         CALL timing_start('tvd_finaltr')
288#if defined key_z_first
289         DO jj = 2, jpjm1
290            DO ji = 2, jpim1
291               DO jk = 1, jpkm1
292#else
293         DO jk = 1, jpkm1
294            DO jj = 2, jpjm1
295               DO ji = fs_2, fs_jpim1   ! vector opt. 
296#endif
297                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
298                  ! total advective trends
299                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
300                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
301                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) )
302                  ! add them to the general tracer trends
303                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
304               END DO
305            END DO
306         END DO
307
308         !                                 ! trend diagnostics (contribution of upstream fluxes)
309         IF( l_trd )  THEN
310            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
311            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
312            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
313           
314            CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )   
315            CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
316            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
317         END IF
318         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
319         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
320           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
321           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
322         ENDIF
323         !
324         CALL timing_stop('tvd_finaltr','section')
325
326      END DO
327      !
328      IF( l_trd )  THEN
329        DEALLOCATE( ztrdx )     ;     DEALLOCATE( ztrdy )     ;      DEALLOCATE( ztrdz ) 
330      END IF
331      !
332      IF( wrk_not_released(3, 12,13) )   CALL ctl_stop('tra_adv_tvd: failed to release workspace arrays')
333      !
334      CALL timing_stop('tra_adv_tvd','section')
335
336      !! * Reset control of array index permutation
337!FTRANS CLEAR
338#  include "oce_ftrans.h90"
339#  include "dom_oce_ftrans.h90"
340#  include "trc_oce_ftrans.h90"
341
342   END SUBROUTINE tra_adv_tvd
343
344
345   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
346      !!---------------------------------------------------------------------
347      !!                    ***  ROUTINE nonosc  ***
348      !!     
349      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
350      !!       scheme and the before field by a nonoscillatory algorithm
351      !!
352      !! **  Method  :   ... ???
353      !!       warning : pbef and paft must be masked, but the boundaries
354      !!       conditions on the fluxes are not necessary zalezak (1979)
355      !!       drange (1995) multi-dimensional forward-in-time and upstream-
356      !!       in-space based differencing for fluid
357      !!----------------------------------------------------------------------
358      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
359      USE wrk_nemo, ONLY:   zbetup => wrk_3d_8  , zbetdo => wrk_3d_9    ! 3D workspace
360      USE wrk_nemo, ONLY:   zbup   => wrk_3d_10 , zbdo   => wrk_3d_11   !  -     -
361      USE timing,   ONLY:   timing_start, timing_stop
362
363      !! DCSE_NEMO: need additional directives for renamed module variables
364!FTRANS zbetup zbetdo zbup zbdo :I :I :z
365
366      !
367      REAL(wp), DIMENSION(jpk)         , INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
368
369      !! DCSE_NEMO: This style defeats ftrans
370!     REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
371!     REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions
372
373!FTRANS pbef paft :I :I :z
374!FTRANS paa pbb pcc :I :I :z
375      REAL(wp), INTENT(in   ) ::   pbef(jpi,jpj,jpk), paft(jpi,jpj,jpk)     ! before & after field
376      REAL(wp), INTENT(inout) ::   paa(jpi,jpj,jpk)                         ! monotonic fluxes in the 1st direction
377      REAL(wp), INTENT(inout) ::   pbb(jpi,jpj,jpk)                         ! monotonic fluxes in the 2nd direction
378      REAL(wp), INTENT(inout) ::   pcc(jpi,jpj,jpk)                         ! monotonic fluxes in the 3rd direction
379      !
380      INTEGER ::   ji, jj, jk   ! dummy loop indices
381      INTEGER ::   ikm1         ! local integer
382      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
383      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
384      !!----------------------------------------------------------------------
385
386      IF( wrk_in_use(3, 8,9,10,11) ) THEN
387         CALL ctl_stop('nonosc: requested workspace array unavailable')   ;   RETURN
388      ENDIF
389
390      zbig  = 1.e+40_wp
391      zrtrn = 1.e-15_wp
392      zbetup(:,:,jpk) = 0._wp   ;   zbetdo(:,:,jpk) = 0._wp
393
394
395      ! Search local extrema
396      ! --------------------
397      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
398      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   &
399         &        paft * tmask - zbig * ( 1.e0 - tmask )  )
400      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   &
401         &        paft * tmask + zbig * ( 1.e0 - tmask )  )
402
403#if defined key_z_first
404      DO jj = 2, jpjm1
405         DO ji = 2, jpim1
406            DO jk = 1, jpkm1
407               ikm1 = MAX(jk-1,1)
408               z2dtt = p2dt(jk)
409#else
410      DO jk = 1, jpkm1
411         ikm1 = MAX(jk-1,1)
412         z2dtt = p2dt(jk)
413         DO jj = 2, jpjm1
414            DO ji = fs_2, fs_jpim1   ! vector opt.
415#endif
416
417               ! search maximum in neighbourhood
418               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
419                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
420                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   &
421                  &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  )
422
423               ! search minimum in neighbourhood
424               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
425                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
426                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   &
427                  &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  )
428
429               ! positive part of the flux
430               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
431                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
432                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
433
434               ! negative part of the flux
435               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
436                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
437                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
438
439               ! up & down beta terms
440               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
441               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
442               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
443            END DO
444         END DO
445      END DO
446      CALL timing_start('tvd_lbc')
447      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
448      CALL timing_stop('tvd_lbc','section')
449
450
451
452      ! 3. monotonic flux in the i & j direction (paa & pbb)
453      ! ----------------------------------------
454#if defined key_z_first
455      DO jj = 2, jpjm1
456         DO ji = 2, jpim1
457            DO jk = 1, jpkm1
458#else
459      DO jk = 1, jpkm1
460         DO jj = 2, jpjm1
461            DO ji = fs_2, fs_jpim1   ! vector opt.
462#endif
463               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
464               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
465               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
466               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
467
468               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
469               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
470               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
471               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
472
473      ! monotonic flux in the k direction, i.e. pcc
474      ! -------------------------------------------
475               za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
476               zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
477               zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
478               pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb )
479            END DO
480         END DO
481      END DO
482      CALL timing_start('tvd_lbc')
483      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
484      CALL timing_stop('tvd_lbc','section')
485
486      !
487      IF( wrk_not_released(3, 8,9,10,11) )   CALL ctl_stop('nonosc: failed to release workspace arrays')
488      !
489   END SUBROUTINE nonosc
490
491   !!======================================================================
492END MODULE traadv_tvd
Note: See TracBrowser for help on using the repository browser.