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

Last change on this file since 408 was 408, checked in by opalod, 18 years ago

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