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.
traadv_tvd.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 258

Last change on this file since 258 was 258, checked in by opalod, 19 years ago

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