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 NEMO/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/traadv_ubs.F90 @ 14903

Last change on this file since 14903 was 14834, checked in by hadcv, 3 years ago

#2600: Merge in dev_r14273_HPC-02_Daley_Tiling
  • Property svn:keywords set to Id
File size: 19.9 KB
Line 
1MODULE traadv_ubs
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_ubs  ***
4   !! Ocean active tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  1.0  !  2006-08  (L. Debreu, R. Benshila)  Original code
7   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
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 oce            ! ocean dynamics and active tracers
15   USE dom_oce        ! ocean space and time domain
16   USE trc_oce        ! share passive tracers/Ocean variables
17   USE trd_oce        ! trends: ocean variables
18   USE traadv_fct      ! acces to routine interp_4th_cpt
19   USE trdtra         ! trends manager: tracers
20   USE diaptr         ! poleward transport diagnostics
21   USE diaar5         ! AR5 diagnostics
22   !
23   USE iom            ! I/O library
24   USE in_out_manager ! I/O manager
25   USE lib_mpp        ! massively parallel library
26   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
27   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
28#if defined key_loop_fusion
29   USE traadv_ubs_lf  ! UBS      scheme            (tra_adv_ubs  routine - loop fusion version)
30#endif
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   tra_adv_ubs   ! routine called by traadv module
36
37   LOGICAL :: l_trd   ! flag to compute trends
38   LOGICAL :: l_ptr   ! flag to compute poleward transport
39   LOGICAL :: l_hst   ! flag to compute heat transport
40
41
42   !! * Substitutions
43#  include "do_loop_substitute.h90"
44#  include "domzgr_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
47   !! $Id$
48   !! Software governed by the CeCILL license (see ./LICENSE)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pU, pV, pW,          &
53      &                    Kbb, Kmm, pt, kjpt, Krhs, kn_ubs_v )
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_adv_ubs  ***
56      !!
57      !! ** Purpose :   Compute the now trend due to the advection of tracers
58      !!      and add it to the general trend of passive tracer equations.
59      !!
60      !! ** Method  :   The 3rd order Upstream Biased Scheme (UBS) is based on an
61      !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
62      !!      It is only used in the horizontal direction.
63      !!      For example the i-component of the advective fluxes are given by :
64      !!                !  e2u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
65      !!          ztu = !  or
66      !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
67      !!      where zltu is the second derivative of the before temperature field:
68      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
69      !!        This results in a dissipatively dominant (i.e. hyper-diffusive)
70      !!      truncation error. The overall performance of the advection scheme
71      !!      is similar to that reported in (Farrow and Stevens, 1995).
72      !!        For stability reasons, the first term of the fluxes which corresponds
73      !!      to a second order centered scheme is evaluated using the now velocity
74      !!      (centered in time) while the second term which is the diffusive part
75      !!      of the scheme, is evaluated using the before velocity (forward in time).
76      !!      Note that UBS is not positive. Do not use it on passive tracers.
77      !!                On the vertical, the advection is evaluated using a FCT scheme,
78      !!      as the UBS have been found to be too diffusive.
79      !!                kn_ubs_v argument controles whether the FCT is based on
80      !!      a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact
81      !!      scheme (kn_ubs_v=4).
82      !!
83      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
84      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
85      !!             - poleward advective heat and salt transport (ln_diaptr=T)
86      !!
87      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
88      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731�1741.
89      !!----------------------------------------------------------------------
90      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
91      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
92      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
93      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
94      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
95      INTEGER                                  , INTENT(in   ) ::   kn_ubs_v        ! number of tracers
96      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
97      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
98      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
99      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
100      !
101      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
102      REAL(wp) ::   ztra, zbtr, zcoef                       ! local scalars
103      REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
104      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
105      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   ztu, ztv, zltu, zltv, zti, ztw     ! 3D workspace
106      !!----------------------------------------------------------------------
107      !
108#if defined key_loop_fusion
109      CALL tra_adv_ubs_lf    ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs, kn_ubs_v )
110#else
111      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
112         IF( kt == kit000 )  THEN
113            IF(lwp) WRITE(numout,*)
114            IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
115            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
116         ENDIF
117         !
118         l_trd = .FALSE.
119         l_hst = .FALSE.
120         l_ptr = .FALSE.
121         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
122         IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.
123         IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
124            &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
125      ENDIF
126      !
127      ztw (:,:, 1 ) = 0._wp      ! surface & bottom value : set to zero for all tracers
128      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp
129      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp
130      !                                                          ! ===========
131      DO jn = 1, kjpt                                            ! tracer loop
132         !                                                       ! ===========
133         !
134         DO jk = 1, jpkm1                !==  horizontal laplacian of before tracer ==!
135            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )                   ! First derivative (masked gradient)
136               zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
137               zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
138               ztu(ji,jj,jk) = zeeu * ( pt(ji+1,jj  ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
139               ztv(ji,jj,jk) = zeev * ( pt(ji  ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
140            END_2D
141            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                   ! Second derivative (divergence)
142               zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) )
143               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
144               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
145            END_2D
146            !
147         END DO
148         IF (nn_hls==1) CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp, ld4only= .TRUE. )   ! Lateral boundary cond. (unchanged sgn)
149         !
150         DO_3D( 1, 0, 1, 0, 1, jpkm1 )   !==  Horizontal advective fluxes  ==!     (UBS)
151            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )        ! upstream transport (x2)
152            zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
153            zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
154            zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
155            !                                                  ! 2nd order centered advective fluxes (x2)
156            zcenut = pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
157            zcenvt = pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
158            !                                                  ! UBS advective fluxes
159            ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
160            ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
161         END_3D
162         !
163         DO_3D( 1, 1, 1, 1, 1, jpk )
164            zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs)      ! store the initial trends before its update
165         END_3D
166         !
167         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
168            DO_2D( 0, 0, 0, 0 )
169               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        &
170                  &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
171                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) &
172                  &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
173            END_2D
174            !
175         END DO
176         !
177         DO_3D( 1, 1, 1, 1, 1, jpk )
178            zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltu(ji,jj,jk)  ! Horizontal advective trend used in vertical 2nd order FCT case
179         END_3D                                                     ! and/or in trend diagnostic (l_trd=T)
180         !
181         IF( l_trd ) THEN                  ! trend diagnostics
182             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztu, pU, pt(:,:,:,jn,Kmm) )
183             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztv, pV, pt(:,:,:,jn,Kmm) )
184         END IF
185         !
186         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
187         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
188         !                                !  heati/salt transport
189         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
190         !
191         !
192         !                       !== vertical advective trend  ==!
193         !
194         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme
195         !
196         CASE(  2  )                   ! 2nd order FCT
197            !
198            IF( l_trd ) THEN
199               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
200                  zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag.
201               END_3D
202            ENDIF
203            !
204            !                               !*  upstream advection with initial mass fluxes & intermediate update  ==!
205            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
206               zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
207               zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
208               ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb)  ) * wmask(ji,jj,jk)
209            END_3D
210            IF( ln_linssh ) THEN                ! top ocean value (only in linear free surface as ztw has been w-masked)
211               IF( ln_isfcav ) THEN                   ! top of the ice-shelf cavities and at the ocean surface
212                  DO_2D( 1, 1, 1, 1 )
213                     ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface
214                  END_2D
215               ELSE                                   ! no cavities: only at the ocean surface
216                  DO_2D( 1, 1, 1, 1 )
217                     ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb)
218                  END_2D
219               ENDIF
220            ENDIF
221            !
222            DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* trend and after field with monotonic scheme
223               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
224                  &     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
225               pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak
226               zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
227            END_3D
228            !
229            !                          !*  anti-diffusive flux : high order minus low order
230            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
231               ztw(ji,jj,jk) = (   0.5_wp * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
232                  &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
233            END_3D
234            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0
235            IF( ln_linssh )   ztw(:,:, 1 ) = 0._wp       ! only ocean surface as interior zwz values have been w-masked
236            !
237            CALL nonosc_z( Kmm, pt(:,:,:,jn,Kbb), ztw, zti, p2dt )      !  monotonicity algorithm
238            !
239         CASE(  4  )                               ! 4th order COMPACT
240            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )         ! 4th order compact interpolation of T at w-point
241            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
242               ztw(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
243            END_3D
244            IF( ln_linssh ) THEN
245               DO_2D( 1, 1, 1, 1 )
246                  ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)     !!gm ISF & 4th COMPACT doesn't work
247               END_2D
248            ENDIF
249            !
250         END SELECT
251         !
252         DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !  final trend with corrected fluxes
253            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
254               &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
255         END_3D
256         !
257         IF( l_trd )  THEN               ! vertical advective trend diagnostics
258            DO_3D( 0, 0, 0, 0, 1, jpkm1 )                 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
259               zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
260                  &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   &
261                  &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
262            END_3D
263            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zltv )
264         ENDIF
265         !
266      END DO
267      !
268#endif
269   END SUBROUTINE tra_adv_ubs
270
271
272   SUBROUTINE nonosc_z( Kmm, pbef, pcc, paft, p2dt )
273      !!---------------------------------------------------------------------
274      !!                    ***  ROUTINE nonosc_z  ***
275      !!
276      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
277      !!       scheme and the before field by a nonoscillatory algorithm
278      !!
279      !! **  Method  :   ... ???
280      !!       warning : pbef and paft must be masked, but the boundaries
281      !!       conditions on the fluxes are not necessary zalezak (1979)
282      !!       drange (1995) multi-dimensional forward-in-time and upstream-
283      !!       in-space based differencing for fluid
284      !!----------------------------------------------------------------------
285      INTEGER , INTENT(in   )                         ::   Kmm    ! time level index
286      REAL(wp), INTENT(in   )                         ::   p2dt   ! tracer time-step
287      REAL(wp),                DIMENSION(jpi,jpj,jpk) ::   pbef   ! before field
288      REAL(wp), INTENT(inout), DIMENSION(A2D(nn_hls)    ,jpk) ::   paft   ! after field
289      REAL(wp), INTENT(inout), DIMENSION(A2D(nn_hls)    ,jpk) ::   pcc    ! monotonic flux in the k direction
290      !
291      INTEGER  ::   ji, jj, jk   ! dummy loop indices
292      INTEGER  ::   ikm1         ! local integer
293      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn   ! local scalars
294      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zbetup, zbetdo         ! 3D workspace
295      !!----------------------------------------------------------------------
296      !
297      zbig  = 1.e+38_wp
298      zrtrn = 1.e-15_wp
299      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
300      !
301      ! Search local extrema
302      ! --------------------
303      !                    ! large negative value (-zbig) inside land
304      DO_3D( 0, 0, 0, 0, 1, jpk )
305         pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1.e0 - tmask(ji,jj,jk) )
306         paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1.e0 - tmask(ji,jj,jk) )
307      END_3D
308      !
309      DO jk = 1, jpkm1     ! search maximum in neighbourhood
310         ikm1 = MAX(jk-1,1)
311         DO_2D( 0, 0, 0, 0 )
312            zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
313               &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
314               &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
315         END_2D
316      END DO
317      !                    ! large positive value (+zbig) inside land
318      DO_3D( 0, 0, 0, 0, 1, jpk )
319         pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1.e0 - tmask(ji,jj,jk) )
320         paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1.e0 - tmask(ji,jj,jk) )
321      END_3D
322      !
323      DO jk = 1, jpkm1     ! search minimum in neighbourhood
324         ikm1 = MAX(jk-1,1)
325         DO_2D( 0, 0, 0, 0 )
326            zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
327               &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
328               &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
329         END_2D
330      END DO
331      !                    ! restore masked values to zero
332      DO_3D( 0, 0, 0, 0, 1, jpk )
333         pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk)
334         paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk)
335      END_3D
336      !
337      ! Positive and negative part of fluxes and beta terms
338      ! ---------------------------------------------------
339      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
340         ! positive & negative part of the flux
341         zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
342         zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
343         ! up & down beta terms
344         zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
345         zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
346         zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
347      END_3D
348      !
349      ! monotonic flux in the k direction, i.e. pcc
350      ! -------------------------------------------
351      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
352         za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
353         zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
354         zc = 0.5 * ( 1.e0 + SIGN( 1.0_wp, pcc(ji,jj,jk) ) )
355         pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
356      END_3D
357      !
358   END SUBROUTINE nonosc_z
359
360   !!======================================================================
361END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.