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/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90 @ 10802

Last change on this file since 10802 was 10802, checked in by davestorkey, 5 years ago

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