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_ubs.F90 in branches/dev_001_GM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 786

Last change on this file since 786 was 786, checked in by gm, 16 years ago

dev_001_GM - merge TRC-TRA on OPA only, trabbl & zpshde not done and trdmld not OK - compilation OK
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 20.6 KB
Line 
1MODULE traadv_ubs
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_ubs  ***
4   !! Ocean tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  1.0  !  06-08  (L. Debreu, R. Benshila)  Original code
7   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   tra_adv_ubs : update the tracer trend with the horizontal
12   !!                 advection trends using a third order biaised scheme 
13   !!----------------------------------------------------------------------
14   USE dom_oce         ! ocean space and time domain
15   USE trdmod
16   USE trdmod_oce
17   USE lib_mpp
18   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
19   USE in_out_manager  ! I/O manager
20   USE diaptr          ! poleward transport diagnostics
21   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
22   USE prtctl
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   tra_adv_ubs   ! routine called by traadv module
28
29   REAL(wp), DIMENSION(jpi,jpj) ::   e1e2tr   ! = 1/(e1t * e2t)
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA & TRP 2.4 , LOCEAN-IPSL (2008)
36   !! $Id:$
37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE tra_adv_ubs( kt, cdtype, ktra, pun, pvn, pwn,   &
43      &                                      ptb, ptn, pta )
44      !!----------------------------------------------------------------------
45      !!                  ***  ROUTINE tra_adv_ubs  ***
46      !!                 
47      !! ** Purpose :   Compute the now trend due to the advection of tracers
48      !!      and add it to the general trend of passive tracer equations.
49      !!
50      !! ** Method  :   The upstream biased third (UBS) is order scheme based
51      !!      on an upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
52      !!      It is only used in the horizontal direction.
53      !!      For example the i-component of the advective fluxes are given by :
54      !!                !  e1u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
55      !!          zwx = !  or
56      !!                !  e1u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
57      !!      where zltu is the second derivative of the before temperature field:
58      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
59      !!      This results in a dissipatively dominant (i.e. hyper-diffusive)
60      !!      truncation error. The overall performance of the advection scheme
61      !!      is similar to that reported in (Farrow and Stevens, 1995).
62      !!      For stability reasons, the first term of the fluxes which corresponds
63      !!      to a second order centered scheme is evaluated using the now velocity
64      !!      (centered in time) while the second term which is the diffusive part
65      !!      of the scheme, is evaluated using the before velocity (forward in time).
66      !!      Note that UBS is not positive. Do not use it on passive tracers.
67      !!      On the vertical, the advection is evaluated using a TVD scheme, as
68      !!      the UBS have been found to be too diffusive.
69      !!
70      !! ** Action : - update (ta,sa) with the now advective tracer trends
71      !!
72      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
73      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731–1741.
74      !!----------------------------------------------------------------------
75      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index
76      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator)
77      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index
78      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! 3 ocean velocity components
79      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn        ! before and now tracer fields
80      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend
81      !!
82      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
83      REAL(wp) ::   zta, zbtr, zcoef                  ! temporary scalars
84      REAL(wp) ::   zfui, zfp_ui, zfm_ui, zcenut      !    "         "
85      REAL(wp) ::   zfvj, zfp_vj, zfm_vj, zcenvt      !    "         "
86      REAL(wp) ::   z_hdivn                           !    "         "
87      REAL(wp), DIMENSION(jpi,jpj)     ::   zeeu, zeev     ! temporary 2D workspace
88      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy              ! temporary 3D workspace
89      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu , ztv , zltu , zltv, ztrdt   !    "              "
90      !!----------------------------------------------------------------------
91
92      zltu(:,:,:) = 0.e0
93      zltv(:,:,:) = 0.e0
94
95      IF( kt == nit000 ) THEN
96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme'
98         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
99         !
100         e1e2tr(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
101      ENDIF
102
103      ! store pta trends
104      ztrdt(:,:,:) = pta(:,:,:)
105
106      zcoef = 1./6.
107      !                                                ! ===============
108      DO jk = 1, jpkm1                                 ! Horizontal slab
109         !                                             ! ===============
110
111         !  Initialization of metric arrays (for z- or s-coordinates)
112         DO jj = 1, jpjm1
113            DO ji = 1, fs_jpim1   ! vector opt.
114#if defined key_zco
115               ! z-coordinates, no vertical scale factors
116               zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
117               zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
118#else
119               ! s-coordinates, vertical scale factor are used
120               zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
121               zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
122#endif
123            END DO
124         END DO
125
126         !  Laplacian
127         ! First derivative (gradient)
128         DO jj = 1, jpjm1
129            DO ji = 1, fs_jpim1   ! vector opt.
130               ztu(ji,jj,jk) = zeeu(ji,jj) * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) )
131               ztv(ji,jj,jk) = zeev(ji,jj) * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) )
132            END DO
133         END DO
134         ! Second derivative (divergence)
135         DO jj = 2, jpjm1
136            DO ji = fs_2, fs_jpim1   ! vector opt.
137#if ! defined key_zco
138               zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
139#endif         
140               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
141               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
142            END DO
143         END DO
144         !                                             ! =================
145      END DO                                           !    End of slab
146      !                                                ! =================
147
148      ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn)
149      CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )
150
151      !                                                ! ===============
152      DO jk = 1, jpkm1                                 ! Horizontal slab
153         !                                             ! ===============
154         !  Horizontal advective fluxes
155         DO jj = 1, jpjm1
156            DO ji = 1, fs_jpim1   ! vector opt.
157               ! volume fluxes * 1/2
158#if defined key_zco
159               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
160               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
161#else
162               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
163               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
164#endif
165               ! upstream scheme
166               zfp_ui = zfui + ABS( zfui )
167               zfp_vj = zfvj + ABS( zfvj )
168               zfm_ui = zfui - ABS( zfui )
169               zfm_vj = zfvj - ABS( zfvj )
170               ! centered scheme
171               zcenut = zfui * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) )
172               zcenvt = zfvj * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) )
173               ! mixed centered / upstream scheme
174               zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk)
175               zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk)
176            END DO
177         END DO
178
179         !  Tracer flux divergence at t-point added to the general trend
180         DO jj = 2, jpjm1
181            DO ji = fs_2, fs_jpim1   ! vector opt.
182               ! horizontal advective trends
183#if defined key_zco
184               zbtr = e1e2tr(ji,jj)
185#else
186               zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk)
187#endif
188               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
189                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
190               ! add it to the general tracer trends
191               pta(ji,jj,jk) = pta(ji,jj,jk) + zta
192            END DO
193         END DO
194         !                                             ! ===============
195      END DO                                           !   End of slab
196      !                                                ! ===============
197
198      ! Horizontal trend used in tra_adv_ztvd subroutine
199      zltu(:,:,:) = pta(:,:,:) - ztrdt(:,:,:) 
200
201      !  Save the horizontal advective trends for diagnostic
202      ! -----------------------------------------------------
203      IF( l_trdtra ) THEN
204         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn )
205         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn )
206      ENDIF
207
208      ! "Poleward" heat or salt transport
209      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
210         IF( lk_zco ) THEN
211            DO jk = 1, jpkm1
212               DO jj = 2, jpjm1
213                  DO ji = fs_2, fs_jpim1   ! vector opt.
214                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
215                  END DO
216               END DO
217            END DO
218         ENDIF
219         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) )
220         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) )
221      ENDIF
222
223      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - had: ', mask1=tmask, clinfo3=cdtype )
224
225
226      ! II. Vertical advection
227      ! ----------------------
228      IF( l_trdtra )   ztrdt(:,:,:) = pta(:,:,:)          ! Save ta and sa trends
229   
230      ! TVD scheme the vertical direction 
231      CALL tra_adv_ztvd( kt, pwn, zltu, ptb, ptn, pta )
232
233      IF( l_trdtra )   THEN         !  vertical advective trend diagnostics
234         DO jk = 1, jpkm1
235            DO jj = 2, jpjm1
236               DO ji = fs_2, fs_jpim1   ! vector opt.
237                  z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) / fse3t(ji,jj,jk)
238                  ztrdt(ji,jj,jk) = pta(ji,jj,jk) - ztrdt(ji,jj,jk+1)  +  ptn(ji,jj,jk) * z_hdivn 
239               END DO
240            END DO
241         END DO
242         CALL trd_tra( kt, ktra, jpt_trd_zad, cdtype, ptrd3d=ztrdt )
243         !
244      ENDIF
245     
246      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - zad: ', mask1=tmask, clinfo3=cdtype )
247      !
248   END SUBROUTINE tra_adv_ubs
249
250
251   SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, ptb, ptn, pta )
252      !!----------------------------------------------------------------------
253      !!                  ***  ROUTINE tra_adv_ztvd  ***
254      !!
255      !! **  Purpose :   Compute the now trend due to total advection of
256      !!       tracers and add it to the general trend of tracer equations
257      !!
258      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with
259      !!       corrected flux (monotonic correction)
260      !!       note: - this advection scheme needs a leap-frog time scheme
261      !!
262      !! ** Action : - update (ta,sa) with the now advective tracer trends
263      !!             - save the trends in (ztrdt,ztrds) ('key_trdtra')
264      !!----------------------------------------------------------------------
265      INTEGER , INTENT(in   )                         ::   kt      ! ocean time-step
266      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn     ! verical effective velocity
267      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   zttrd   ! lateral advective trends on T & S
268      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn        ! before and now tracer fields
269      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend
270      !!
271      INTEGER  ::   ji, jj, jk              ! dummy loop indices
272      REAL(wp) ::   z2dtt, zbtr, zew, z2    ! temporary scalar 
273      REAL(wp) ::   ztak, zfp_wk, zfm_wk            !    "         "
274      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztw   ! temporary 3D workspace
275      !!----------------------------------------------------------------------
276
277      IF( kt == nit000 .AND. lwp ) THEN
278         WRITE(numout,*)
279         WRITE(numout,*) 'tra_adv_ztvd : vertical TVD advection scheme'
280         WRITE(numout,*) '~~~~~~~~~~~~'
281      ENDIF
282
283      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
284      ELSE                                        ;    z2 = 2.
285      ENDIF
286
287      !  Bottom value : flux set to zero
288      ! --------------
289      ztw(:,:,jpk) = 0.e0   ;   zti  (:,:,:) = 0.e0
290
291
292      !  upstream advection with initial mass fluxes & intermediate update
293      ! -------------------------------------------------------------------
294      ! Surface value
295      IF( lk_dynspg_rl .OR. lk_vvl ) THEN                           ! rigid lid : flux set to zero
296         ztw(:,:,1) = 0.e0
297      ELSE                                                          ! free surface
298         ztw(:,:,1) = e1t(:,:) * e2t(:,:) * pwn(:,:,1) * ptb(:,:,1)
299      ENDIF
300
301      ! Interior value
302      DO jk = 2, jpkm1
303         DO jj = 1, jpj
304            DO ji = 1, jpi
305               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
306               zfp_wk = zew + ABS( zew )
307               zfm_wk = zew - ABS( zew )
308               ztw(ji,jj,jk) = zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1)
309            END DO
310         END DO
311      END DO
312
313      ! update and guess with monotonic sheme
314      DO jk = 1, jpkm1
315         z2dtt = z2 * rdttra(jk)
316         DO jj = 2, jpjm1
317            DO ji = fs_2, fs_jpim1   ! vector opt.
318               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
319               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
320               pta(ji,jj,jk) =  pta(ji,jj,jk) + ztak
321               zti (ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * ( ztak + zttrd(ji,jj,jk) ) ) * tmask(ji,jj,jk)
322            END DO
323         END DO
324      END DO
325
326      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
327      CALL lbc_lnk( zti, 'T', 1. )
328
329
330      !  antidiffusive flux : high order minus low order
331      ! -------------------------------------------------     
332      ztw(:,:,1) = 0.e0       ! Surface value
333
334      DO jk = 2, jpkm1        ! Interior value
335         DO jj = 1, jpj
336            DO ji = 1, jpi
337               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk)
338               ztw(ji,jj,jk) = zew * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
339            END DO
340         END DO
341      END DO
342
343      !  monotonicity algorithm
344      ! ------------------------
345      CALL nonosc_z( ptb, ztw, zti, z2 )
346
347
348      !  final trend with corrected fluxes
349      ! -----------------------------------
350      DO jk = 1, jpkm1
351         DO jj = 2, jpjm1
352            DO ji = fs_2, fs_jpim1   ! vector opt. 
353               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
354               ! k- vertical advective trends
355               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
356               ! add them to the general tracer trends
357               pta(ji,jj,jk) = pta(ji,jj,jk) + ztak
358            END DO
359         END DO
360      END DO
361      !
362   END SUBROUTINE tra_adv_ztvd
363
364
365   SUBROUTINE nonosc_z( pbef, pcc, paft, prdt )
366      !!---------------------------------------------------------------------
367      !!                    ***  ROUTINE nonosc_z  ***
368      !!     
369      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
370      !!       scheme and the before field by a nonoscillatory algorithm
371      !!
372      !! **  Method  :   ... ???
373      !!       warning : pbef and paft must be masked, but the boundaries
374      !!       conditions on the fluxes are not necessary zalezak (1979)
375      !!       drange (1995) multi-dimensional forward-in-time and upstream-
376      !!       in-space based differencing for fluid
377      !!----------------------------------------------------------------------
378      REAL(wp), INTENT(in   )                          ::   prdt   ! ???
379      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
380      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
381      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
382      !!
383      INTEGER  ::   ji, jj, jk               ! dummy loop indices
384      INTEGER  ::   ikm1
385      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
386      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
387      !!----------------------------------------------------------------------
388
389      zbig = 1.e+40
390      zrtrn = 1.e-15
391      zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0
392
393      ! Search local extrema
394      ! --------------------
395      ! large negative value (-zbig) inside land
396      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
397      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
398      ! search maximum in neighbourhood
399      DO jk = 1, jpkm1
400         ikm1 = MAX(jk-1,1)
401         DO jj = 2, jpjm1
402            DO ji = fs_2, fs_jpim1   ! vector opt.
403               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
404                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
405                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
406            END DO
407         END DO
408      END DO
409      ! large positive value (+zbig) inside land
410      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
411      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
412      ! search minimum in neighbourhood
413      DO jk = 1, jpkm1
414         ikm1 = MAX(jk-1,1)
415         DO jj = 2, jpjm1
416            DO ji = fs_2, fs_jpim1   ! vector opt.
417               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
418                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
419                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
420            END DO
421         END DO
422      END DO
423
424      ! restore masked values to zero
425      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
426      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
427 
428
429      ! 2. Positive and negative part of fluxes and beta terms
430      ! ------------------------------------------------------
431
432      DO jk = 1, jpkm1
433         z2dtt = prdt * rdttra(jk)
434         DO jj = 2, jpjm1
435            DO ji = fs_2, fs_jpim1   ! vector opt.
436               ! positive & negative part of the flux
437               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
438               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
439               ! up & down beta terms
440               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
441               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
442               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
443            END DO
444         END DO
445      END DO
446
447      ! monotonic flux in the k direction, i.e. pcc
448      ! -------------------------------------------
449      DO jk = 2, jpkm1
450         DO jj = 2, jpjm1
451            DO ji = fs_2, fs_jpim1   ! vector opt.
452
453               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
454               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
455               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
456               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
457            END DO
458         END DO
459      END DO
460      !
461   END SUBROUTINE nonosc_z
462
463   !!======================================================================
464END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.