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

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

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