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

Last change on this file since 1528 was 1528, checked in by rblod, 15 years ago

Suppress rigid-lid option, see ticket #486
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 22.2 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   !! $Id$
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_vvl ) THEN
132         ! variable volume : flux set to zero
133         ztw(:,:,1) = 0.e0
134         zsw(:,:,1) = 0.e0
135      ELSE
136         ! free surface-constant volume
137         DO jj = 1, jpj
138            DO ji = 1, jpi
139               zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1)
140               ztw(ji,jj,1) = zew * tb(ji,jj,1)
141               zsw(ji,jj,1) = zew * sb(ji,jj,1)
142            END DO
143         END DO
144      ENDIF
145
146      ! Interior value
147      DO jk = 2, jpkm1
148         DO jj = 1, jpj
149            DO ji = 1, jpi
150               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
151               zfp_wk = zew + ABS( zew )
152               zfm_wk = zew - ABS( zew )
153               ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1)
154               zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1)
155            END DO
156         END DO
157      END DO
158
159      ! total advective trend
160      DO jk = 1, jpkm1
161         z2dtt = z2 * rdttra(jk)
162         DO jj = 2, jpjm1
163            DO ji = fs_2, fs_jpim1   ! vector opt.
164               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
165               ! total intermediate advective trends
166               ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )   &
167                  &     + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )   &
168                  &     + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr
169               zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  )   & 
170                  &     + zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  )   &
171                  &     + zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1) ) * zbtr
172               ! update and guess with monotonic sheme
173               ta(ji,jj,jk) =  ta(ji,jj,jk) + ztat
174               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsat
175               zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ztat ) * tmask(ji,jj,jk)
176               zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * zsat ) * tmask(ji,jj,jk)
177            END DO
178         END DO
179      END DO
180
181
182      ! Save the intermediate i / j / k advective trends for diagnostics
183      ! -------------------------------------------------------------------
184      ! Warning : We should use zun instead of un in the computations below, but we
185      ! also use hdivn which is computed with un, vn (check ???). So we use un, vn
186      ! for consistency. Results are therefore approximate with key_trabbl_adv.
187
188      IF( l_trdtra ) THEN
189         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
190         !
191         ! T/S ZONAL advection trends
192         DO jk = 1, jpkm1
193            DO jj = 2, jpjm1
194               DO ji = fs_2, fs_jpim1   ! vector opt.
195                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
196                  ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr
197                  ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr
198               END DO
199            END DO
200         END DO
201         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)    ! save the trends
202         !
203         ! T/S MERIDIONAL advection trends
204         DO jk = 1, jpkm1
205            DO jj = 2, jpjm1
206               DO ji = fs_2, fs_jpim1   ! vector opt.
207                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
208                  ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr
209                  ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr
210               END DO
211            END DO
212         END DO
213         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)     ! save the trends
214         !
215         ! T/S VERTICAL advection trends
216         DO jk = 1, jpkm1
217            DO jj = 2, jpjm1
218               DO ji = fs_2, fs_jpim1   ! vector opt.         
219                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
220                  ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
221                  ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr
222               END DO
223            END DO
224         END DO
225         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)     ! save the trends
226         !
227      ENDIF
228
229      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
230      CALL lbc_lnk( zti, 'T', 1. )
231      CALL lbc_lnk( zsi, 'T', 1. )
232
233
234      ! 3. antidiffusive flux : high order minus low order
235      ! --------------------------------------------------
236      ! antidiffusive flux on i and j
237      DO jk = 1, jpkm1
238         DO jj = 1, jpjm1
239            DO ji = 1, fs_jpim1   ! vector opt.
240               zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
241               zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
242               ztu(ji,jj,jk) = zeu * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) ) - ztu(ji,jj,jk)
243               zsu(ji,jj,jk) = zeu * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) ) - zsu(ji,jj,jk)
244               ztv(ji,jj,jk) = zev * ( tn(ji,jj,jk) + tn(ji,jj+1,jk) ) - ztv(ji,jj,jk)
245               zsv(ji,jj,jk) = zev * ( sn(ji,jj,jk) + sn(ji,jj+1,jk) ) - zsv(ji,jj,jk)
246            END DO
247         END DO
248      END DO
249     
250      ! antidiffusive flux on k
251      ! Surface value
252      ztw(:,:,1) = 0.e0
253      zsw(:,:,1) = 0.e0
254
255      ! Interior value
256      DO jk = 2, jpkm1
257         DO jj = 1, jpj
258            DO ji = 1, jpi
259               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
260               ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
261               zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk)
262            END DO
263         END DO
264      END DO
265
266      ! Lateral bondary conditions
267      CALL lbc_lnk( ztu, 'U', -1. )   ;   CALL lbc_lnk( zsu, 'U', -1. )
268      CALL lbc_lnk( ztv, 'V', -1. )   ;   CALL lbc_lnk( zsv, 'V', -1. )
269      CALL lbc_lnk( ztw, 'W',  1. )   ;   CALL lbc_lnk( zsw, 'W',  1. )
270
271      ! 4. monotonicity algorithm
272      ! -------------------------
273      CALL nonosc( tb, ztu, ztv, ztw, zti, z2 )
274      CALL nonosc( sb, zsu, zsv, zsw, zsi, z2 )
275
276
277      ! 5. final trend with corrected fluxes
278      ! ------------------------------------
279      DO jk = 1, jpkm1
280         DO jj = 2, jpjm1
281            DO ji = fs_2, fs_jpim1   ! vector opt. 
282               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
283               ! total advective trends
284               ztat = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )   &
285                  &     + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )   &
286                  &     + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr
287               zsat = - ( zsu(ji,jj,jk) - zsu(ji-1,jj  ,jk  )   &
288                  &     + zsv(ji,jj,jk) - zsv(ji  ,jj-1,jk  )   &
289                  &     + zsw(ji,jj,jk) - zsw(ji  ,jj  ,jk+1) ) * zbtr
290               ! add them to the general tracer trends
291               ta(ji,jj,jk) = ta(ji,jj,jk) + ztat
292               sa(ji,jj,jk) = sa(ji,jj,jk) + zsat
293            END DO
294         END DO
295      END DO
296
297
298      ! Save the advective trends for diagnostics
299      ! --------------------------------------------
300
301      IF( l_trdtra ) THEN
302         ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
303         !
304         ! T/S ZONAL advection trends
305         DO jk = 1, jpkm1
306            DO jj = 2, jpjm1
307               DO ji = fs_2, fs_jpim1   ! vector opt.
308                  !-- Compute zonal divergence by splitting hdivn (see divcur.F90)
309                  !   N.B. This computation is not valid along OBCs (if any)
310                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
311                  z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
312                     &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
313                  !-- Compute T/S zonal advection trends
314                  ztrdt(ji,jj,jk) = - ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x
315                  ztrds(ji,jj,jk) = - ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x
316               END DO
317            END DO
318         END DO
319         CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt, cnbpas='bis')   ! <<< ADD TO PREVIOUSLY COMPUTED
320         !
321         ! T/S MERIDIONAL advection trends
322         DO jk = 1, jpkm1
323            DO jj = 2, jpjm1
324               DO ji = fs_2, fs_jpim1   ! vector opt.
325                  !-- Compute merid. divergence by splitting hdivn (see divcur.F90)
326                  !   N.B. This computation is not valid along OBCs (if any)
327                  zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
328                  z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          &
329                     &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
330                  !-- Compute T/S meridional advection trends
331                  ztrdt(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y         
332                  ztrds(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y         
333               END DO
334            END DO
335         END DO
336         CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt, cnbpas='bis')   ! <<< ADD TO PREVIOUSLY COMPUTED
337         !
338         ! T/S VERTICAL advection trends
339         DO jk = 1, jpkm1
340            DO jj = 2, jpjm1
341               DO ji = fs_2, fs_jpim1   ! vector opt.
342                  zbtr1     = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
343#if defined key_zco
344                  zbtr      = zbtr1
345                  z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk)
346                  z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk)
347#else
348                  zbtr      = zbtr1 / fse3t(ji,jj,jk)
349                  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)
350                  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)
351#endif
352                  z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr
353                  zbtr      = zbtr1 / fse3t(ji,jj,jk)
354                  ztrdt(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr - tn(ji,jj,jk) * z_hdivn
355                  ztrds(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr - sn(ji,jj,jk) * z_hdivn
356               END DO
357            END DO
358         END DO
359         CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt, cnbpas='bis')   ! <<< ADD TO PREVIOUSLY COMPUTED
360         !
361      ENDIF
362
363      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' tvd adv  - Ta: ', mask1=tmask,   &
364         &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
365
366      ! "zonal" mean advective heat and salt transport
367      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
368         pht_adv(:) = ptr_vj( ztv(:,:,:) )
369         pst_adv(:) = ptr_vj( zsv(:,:,:) )
370      ENDIF
371      !
372   END SUBROUTINE tra_adv_tvd
373
374
375   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, prdt )
376      !!---------------------------------------------------------------------
377      !!                    ***  ROUTINE nonosc  ***
378      !!     
379      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
380      !!       scheme and the before field by a nonoscillatory algorithm
381      !!
382      !! **  Method  :   ... ???
383      !!       warning : pbef and paft must be masked, but the boundaries
384      !!       conditions on the fluxes are not necessary zalezak (1979)
385      !!       drange (1995) multi-dimensional forward-in-time and upstream-
386      !!       in-space based differencing for fluid
387      !!----------------------------------------------------------------------
388      REAL(wp), INTENT( in ) ::   prdt   ! ???
389      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( inout ) ::   &
390         pbef,                            & ! before field
391         paft,                            & ! after field
392         paa,                             & ! monotonic flux in the i direction
393         pbb,                             & ! monotonic flux in the j direction
394         pcc                                ! monotonic flux in the k direction
395      !!
396      INTEGER ::   ji, jj, jk               ! dummy loop indices
397      INTEGER ::   ikm1
398      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
399      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbup, zbdo
400      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
401      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv
402      REAL(wp) ::   zup, zdo
403      !!----------------------------------------------------------------------
404
405      zbig = 1.e+40
406      zrtrn = 1.e-15
407      zbetup(:,:,jpk) = 0.e0   ;   zbetdo(:,:,jpk) = 0.e0
408
409
410      ! Search local extrema
411      ! --------------------
412      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
413      zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ),   &
414         &        paft * tmask - zbig * ( 1.e0 - tmask )  )
415      zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ),   &
416         &        paft * tmask + zbig * ( 1.e0 - tmask )  )
417
418      DO jk = 1, jpkm1
419         ikm1 = MAX(jk-1,1)
420         z2dtt = prdt * rdttra(jk)
421         DO jj = 2, jpjm1
422            DO ji = fs_2, fs_jpim1   ! vector opt.
423
424               ! search maximum in neighbourhood
425               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
426                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
427                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   &
428                  &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  )
429
430               ! search minimum in neighbourhood
431               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
432                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
433                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   &
434                  &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  )
435
436               ! positive part of the flux
437               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
438                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
439                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
440
441               ! negative part of the flux
442               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
443                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
444                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
445
446               ! up & down beta terms
447               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
448               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
449               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
450            END DO
451         END DO
452      END DO
453
454      ! lateral boundary condition on zbetup & zbetdo   (unchanged sign)
455      CALL lbc_lnk( zbetup, 'T', 1. )
456      CALL lbc_lnk( zbetdo, 'T', 1. )
457
458
459      ! 3. monotonic flux in the i & j direction (paa & pbb)
460      ! ----------------------------------------
461      DO jk = 1, jpkm1
462         DO jj = 2, jpjm1
463            DO ji = fs_2, fs_jpim1   ! vector opt.
464               zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
465               zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
466               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
467               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
468
469               zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
470               zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
471               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
472               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
473
474      ! monotonic flux in the k direction, i.e. pcc
475      ! -------------------------------------------
476               za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
477               zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
478               zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
479               pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb )
480            END DO
481         END DO
482      END DO
483
484      ! lateral boundary condition on paa, pbb, pcc
485      CALL lbc_lnk( paa, 'U', -1. )      ! changed sign
486      CALL lbc_lnk( pbb, 'V', -1. )      ! changed sign
487      !
488   END SUBROUTINE nonosc
489
490   !!======================================================================
491END MODULE traadv_tvd
Note: See TracBrowser for help on using the repository browser.