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 branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 7236

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

All changes related to diaptr (basin heat transports and transport components)

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