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/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_ubs.F90 @ 13819

Last change on this file since 13819 was 13819, checked in by hadcv, 4 years ago

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