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

Last change on this file since 785 was 785, checked in by rblod, 16 years ago

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