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 @ 237

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

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