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 trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 7881

Last change on this file since 7881 was 7646, checked in by timgraham, 7 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

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