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

Last change on this file since 10880 was 10880, checked in by davestorkey, 3 years ago

branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps:

  1. Move time indices from dom_oce.F90 to step.F90.
  2. Implement time dimension for passive tracers with independent set of time indices in trcstp.F90.
  3. Update all the traadv and trcadv modules.
  • Property svn:keywords set to Id
File size: 20.6 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, cdtype, p2dt, pu_mm, pv_mm, pww,          &
49      &                    Kbb, Kmm, pt, kjpt, Krhs, 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  ) ,Kmm)   if uu(i,Kmm) >= 0
61      !!          ztu = !  or
62      !!                !  e2u e3u uu ( mi(Tn) - zltu(i+1) ,Kmm)   if uu(i,Kmm) < 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(:,:,:,:,Krhs)  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   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
88      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
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_mm, pv_mm, pww   ! 3 ocean transport components
94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
95      !
96      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
97      REAL(wp) ::   ztra, zbtr, zcoef                       ! local scalars
98      REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
99      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
100      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu, ztv, zltu, zltv, zti, ztw   ! 3D workspace
101      !!----------------------------------------------------------------------
102      !
103      IF( kt == kit000 )  THEN
104         IF(lwp) WRITE(numout,*)
105         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
106         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
107      ENDIF
108      !
109      l_trd = .FALSE.
110      l_hst = .FALSE.
111      l_ptr = .FALSE.
112      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
113      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
114      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
115         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
116      !
117      ztw (:,:, 1 ) = 0._wp      ! surface & bottom value : set to zero for all tracers
118      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp
119      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp
120      !
121      !                                                          ! ===========
122      DO jn = 1, kjpt                                            ! tracer loop
123         !                                                       ! ===========
124         !                                             
125         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==!
126            DO jj = 1, jpjm1              ! First derivative (masked gradient)
127               DO ji = 1, fs_jpim1   ! vector opt.
128                  zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
129                  zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
130                  ztu(ji,jj,jk) = zeeu * ( pt(ji+1,jj  ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
131                  ztv(ji,jj,jk) = zeev * ( pt(ji  ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
132               END DO
133            END DO
134            DO jj = 2, jpjm1              ! Second derivative (divergence)
135               DO ji = fs_2, fs_jpim1   ! vector opt.
136                  zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) )
137                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
138                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
139               END DO
140            END DO
141            !                                   
142         END DO         
143         CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1. )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
144         !   
145         DO jk = 1, jpkm1        !==  Horizontal advective fluxes  ==!     (UBS)
146            DO jj = 1, jpjm1
147               DO ji = 1, fs_jpim1   ! vector opt.
148                  zfp_ui = pu_mm(ji,jj,jk) + ABS( pu_mm(ji,jj,jk) )      ! upstream transport (x2)
149                  zfm_ui = pu_mm(ji,jj,jk) - ABS( pu_mm(ji,jj,jk) )
150                  zfp_vj = pv_mm(ji,jj,jk) + ABS( pv_mm(ji,jj,jk) )
151                  zfm_vj = pv_mm(ji,jj,jk) - ABS( pv_mm(ji,jj,jk) )
152                  !                                                  ! 2nd order centered advective fluxes (x2)
153                  zcenut = pu_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
154                  zcenvt = pv_mm(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
155                  !                                                  ! UBS advective fluxes
156                  ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
157                  ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
158               END DO
159            END DO
160         END DO         
161         !
162         zltu(:,:,:) = pt(:,:,:,jn,Krhs)      ! store the initial trends before its update
163         !
164         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
165            DO jj = 2, jpjm1
166               DO ji = fs_2, fs_jpim1   ! vector opt.
167                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        &
168                     &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
169                     &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
170               END DO
171            END DO
172            !                                             
173         END DO
174         !
175         zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
176         !                                            ! and/or in trend diagnostic (l_trd=T)
177         !               
178         IF( l_trd ) THEN                  ! trend diagnostics
179             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pu_mm, pt(:,:,:,jn,Kmm) )
180             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pv_mm, pt(:,:,:,jn,Kmm) )
181         END IF
182         !     
183         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
184         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
185         !                                !  heati/salt transport
186         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
187         !
188         !
189         !                       !== vertical advective trend  ==!
190         !
191         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme
192         !
193         CASE(  2  )                   ! 2nd order FCT
194            !         
195            IF( l_trd )   zltv(:,:,:) = pt(:,:,:,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag.
196            !
197            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
198            DO jk = 2, jpkm1                 ! Interior value (w-masked)
199               DO jj = 1, jpj
200                  DO ji = 1, jpi
201                     zfp_wk = pww(ji,jj,jk) + ABS( pww(ji,jj,jk) )
202                     zfm_wk = pww(ji,jj,jk) - ABS( pww(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 DO
205               END DO
206            END DO
207            IF( ln_linssh ) THEN             ! top ocean value (only in linear free surface as ztw has been w-masked)
208               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
209                  DO jj = 1, jpj
210                     DO ji = 1, jpi
211                        ztw(ji,jj, mikt(ji,jj) ) = pww(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface
212                     END DO
213                  END DO   
214               ELSE                                ! no cavities: only at the ocean surface
215                  ztw(:,:,1) = pww(:,:,1) * pt(:,:,1,jn,Kbb)
216               ENDIF
217            ENDIF
218            !
219            DO jk = 1, jpkm1           !* trend and after field with monotonic scheme
220               DO jj = 2, jpjm1
221                  DO ji = fs_2, fs_jpim1   ! vector opt.
222                     ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
223                     pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak 
224                     zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
225                  END DO
226               END DO
227            END DO
228            CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
229            !
230            !                          !*  anti-diffusive flux : high order minus low order
231            DO jk = 2, jpkm1        ! Interior value  (w-masked)
232               DO jj = 1, jpj
233                  DO ji = 1, jpi
234                     ztw(ji,jj,jk) = (   0.5_wp * pww(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
235                        &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
236                  END DO
237               END DO
238            END DO
239            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0
240            IF( ln_linssh )   ztw(:,:, 1 ) = 0._wp       ! only ocean surface as interior zwz values have been w-masked
241            !
242            CALL nonosc_z( Kmm, pt(:,:,:,jn,Kbb), ztw, zti, p2dt )      !  monotonicity algorithm
243            !
244         CASE(  4  )                               ! 4th order COMPACT
245            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )         ! 4th order compact interpolation of T at w-point
246            DO jk = 2, jpkm1
247               DO jj = 2, jpjm1
248                  DO ji = fs_2, fs_jpim1
249                     ztw(ji,jj,jk) = pww(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
250                  END DO
251               END DO
252            END DO
253            IF( ln_linssh )   ztw(:,:, 1 ) = pww(:,:,1) * pt(:,:,1,jn,Kmm)     !!gm ISF & 4th COMPACT doesn't work
254            !
255         END SELECT
256         !
257         DO jk = 1, jpkm1        !  final trend with corrected fluxes
258            DO jj = 2, jpjm1 
259               DO ji = fs_2, fs_jpim1   ! vector opt.   
260                  pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
261               END DO
262            END DO
263         END DO
264         !
265         IF( l_trd )  THEN       ! vertical advective trend diagnostics
266            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
267               DO jj = 2, jpjm1
268                  DO ji = fs_2, fs_jpim1   ! vector opt.
269                     zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
270                        &           + pt(ji,jj,jk,jn,Kmm) * (  pww(ji,jj,jk) - pww(ji,jj,jk+1)  )   &
271                        &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
272                  END DO
273               END DO
274            END DO
275            CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv )
276         ENDIF
277         !
278      END DO
279      !
280   END SUBROUTINE tra_adv_ubs
281
282
283   SUBROUTINE nonosc_z( Kmm, pbef, pcc, paft, p2dt )
284      !!---------------------------------------------------------------------
285      !!                    ***  ROUTINE nonosc_z  ***
286      !!     
287      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
288      !!       scheme and the before field by a nonoscillatory algorithm
289      !!
290      !! **  Method  :   ... ???
291      !!       warning : pbef and paft must be masked, but the boundaries
292      !!       conditions on the fluxes are not necessary zalezak (1979)
293      !!       drange (1995) multi-dimensional forward-in-time and upstream-
294      !!       in-space based differencing for fluid
295      !!----------------------------------------------------------------------
296      INTEGER , INTENT(in   )                          ::   Kmm    ! time level index
297      REAL(wp), INTENT(in   )                          ::   p2dt   ! tracer time-step
298      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
299      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
300      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
301      !
302      INTEGER  ::   ji, jj, jk   ! dummy loop indices
303      INTEGER  ::   ikm1         ! local integer
304      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn   ! local scalars
305      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zbetup, zbetdo     ! 3D workspace
306      !!----------------------------------------------------------------------
307      !
308      zbig  = 1.e+40_wp
309      zrtrn = 1.e-15_wp
310      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
311      !
312      ! Search local extrema
313      ! --------------------
314      !                    ! large negative value (-zbig) inside land
315      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
316      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
317      !
318      DO jk = 1, jpkm1     ! search maximum in neighbourhood
319         ikm1 = MAX(jk-1,1)
320         DO jj = 2, jpjm1
321            DO ji = fs_2, fs_jpim1   ! vector opt.
322               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
323                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
324                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
325            END DO
326         END DO
327      END DO
328      !                    ! large positive value (+zbig) inside land
329      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
330      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
331      !
332      DO jk = 1, jpkm1     ! search minimum in neighbourhood
333         ikm1 = MAX(jk-1,1)
334         DO jj = 2, jpjm1
335            DO ji = fs_2, fs_jpim1   ! vector opt.
336               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
337                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
338                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
339            END DO
340         END DO
341      END DO
342      !                    ! restore masked values to zero
343      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
344      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
345      !
346      ! Positive and negative part of fluxes and beta terms
347      ! ---------------------------------------------------
348      DO jk = 1, jpkm1
349         DO jj = 2, jpjm1
350            DO ji = fs_2, fs_jpim1   ! vector opt.
351               ! positive & negative part of the flux
352               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
353               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
354               ! up & down beta terms
355               zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
356               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
357               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
358            END DO
359         END DO
360      END DO
361      !
362      ! monotonic flux in the k direction, i.e. pcc
363      ! -------------------------------------------
364      DO jk = 2, jpkm1
365         DO jj = 2, jpjm1
366            DO ji = fs_2, fs_jpim1   ! vector opt.
367               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
368               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
369               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
370               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
371            END DO
372         END DO
373      END DO
374      !
375   END SUBROUTINE nonosc_z
376
377   !!======================================================================
378END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.