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

Last change on this file since 13516 was 13516, checked in by hadcv, 7 months ago

Tiling for tra_adv

  • Property svn:keywords set to Id
File size: 21.2 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   ! TEMP: This change not necessary after trd_tra is tiled
17   USE domain, ONLY : dom_tile
18   USE trc_oce        ! share passive tracers/Ocean variables
19   USE trd_oce        ! trends: ocean variables
20   USE traadv_fct      ! acces to routine interp_4th_cpt
21   USE trdtra         ! trends manager: tracers
22   USE diaptr         ! poleward transport diagnostics
23   USE diaar5         ! AR5 diagnostics
24   !
25   USE iom            ! I/O library
26   USE in_out_manager ! I/O manager
27   USE lib_mpp        ! massively parallel library
28   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
29   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   tra_adv_ubs   ! routine called by traadv module
35
36   LOGICAL :: l_trd   ! flag to compute trends
37   LOGICAL :: l_ptr   ! flag to compute poleward transport
38   LOGICAL :: l_hst   ! flag to compute heat transport
39
40
41   !! * Substitutions
42#  include "do_loop_substitute.h90"
43#  include "domzgr_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
46   !! $Id$
47   !! Software governed by the CeCILL license (see ./LICENSE)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pU, pV, pW,          &
52      &                    Kbb, Kmm, pt, kjpt, Krhs, 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 pt(:,:,:,:,Krhs)  with the now advective tracer trends
83      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
84      !!             - 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   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
91      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
92      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
93      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
94      INTEGER                                  , INTENT(in   ) ::   kn_ubs_v        ! number of tracers
95      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
96      ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
97      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
98      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
99      !
100      ! TEMP: This change not necessary after trd_tra is tiled
101      INTEGER  ::   itile
102      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
103      REAL(wp) ::   ztra, zbtr, zcoef                       ! local scalars
104      REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
105      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
106      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk) ::   ztu, ztv, zltu, zltv, zti, ztw     ! 3D workspace
107      ! TEMP: This change not necessary after trd_tra is tiled
108      REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   ztrdx, ztrdy, ztrdz
109      !!----------------------------------------------------------------------
110      ! TEMP: This change not necessary after trd_tra is tiled
111      itile = ntile
112      !
113      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
114         IF( kt == kit000 )  THEN
115            IF(lwp) WRITE(numout,*)
116            IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
117            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
118         ENDIF
119         !
120         l_trd = .FALSE.
121         l_hst = .FALSE.
122         l_ptr = .FALSE.
123         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
124         IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.
125         IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
126            &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
127
128         ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
129         IF( kt == kit000 .AND. l_trd ) THEN
130            ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) )
131         ENDIF
132      ENDIF
133      !
134      ztw (:,:, 1 ) = 0._wp      ! surface & bottom value : set to zero for all tracers
135      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp
136      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp
137      !
138      !                                                          ! ===========
139      DO jn = 1, kjpt                                            ! tracer loop
140         !                                                       ! ===========
141         !                                             
142         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==!
143            DO_2D( 1, 0, 1, 0 )
144               zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
145               zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
146               ztu(ji,jj,jk) = zeeu * ( pt(ji+1,jj  ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
147               ztv(ji,jj,jk) = zeev * ( pt(ji  ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
148            END_2D
149            DO_2D( 0, 0, 0, 0 )
150               zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) )
151               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
152               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
153            END_2D
154            !                                   
155         END DO         
156         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)
157         !   
158         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
159            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )      ! upstream transport (x2)
160            zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
161            zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
162            zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
163            !                                                  ! 2nd order centered advective fluxes (x2)
164            zcenut = pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
165            zcenvt = pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
166            !                                                  ! UBS advective fluxes
167            ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
168            ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
169         END_3D
170         !
171         DO_3D( 1, 1, 1, 1, 1, jpk )
172            zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs)      ! store the initial trends before its update
173         END_3D
174         !
175         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
176            DO_2D( 0, 0, 0, 0 )
177               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        &
178                  &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
179                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) &
180                  &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
181            END_2D
182            !                                             
183         END DO
184         !
185         DO_3D( 1, 1, 1, 1, 1, jpk )
186            zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltu(ji,jj,jk)  ! Horizontal advective trend used in vertical 2nd order FCT case
187         END_3D                                                     ! and/or in trend diagnostic (l_trd=T)
188         !
189         ! TEMP: These changes not necessary after trd_tra is tiled
190         IF( l_trd ) THEN                  ! trend diagnostics
191            DO_3D( 1, 0, 1, 0, 1, jpkm1 )
192               ztrdx(ji,jj,jk) = ztu(ji,jj,jk)
193               ztrdy(ji,jj,jk) = ztv(ji,jj,jk)
194            END_3D
195
196            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
197               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
198
199               ! TODO: TO BE TILED- trd_tra
200               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) )
201               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) )
202
203               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = itile )     ! Revert to tile domain
204            ENDIF
205         END IF
206         !     
207         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
208         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
209         !                                !  heati/salt transport
210         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
211         !
212         !
213         !                       !== vertical advective trend  ==!
214         !
215         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme
216         !
217         CASE(  2  )                   ! 2nd order FCT
218            !         
219            IF( l_trd ) THEN
220               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
221                  zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag.
222               END_3D
223            ENDIF
224            !
225            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
226            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
227               zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
228               zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
229               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)
230            END_3D
231            IF( ln_linssh ) THEN             ! top ocean value (only in linear free surface as ztw has been w-masked)
232               ! TODO: NOT TESTED- requires isf
233               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
234                  DO_2D( 1, 1, 1, 1 )
235                     ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface
236                  END_2D
237               ELSE                                ! no cavities: only at the ocean surface
238                  DO_2D( 1, 1, 1, 1 )
239                     ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb)
240                  END_2D
241               ENDIF
242            ENDIF
243            !
244            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
245               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
246                  &     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
247               pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak 
248               zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
249            END_3D
250            CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1.0_wp )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
251            !
252            !                          !*  anti-diffusive flux : high order minus low order
253            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
254               ztw(ji,jj,jk) = (   0.5_wp * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
255                  &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
256            END_3D
257            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0
258            IF( ln_linssh )   ztw(:,:, 1 ) = 0._wp       ! only ocean surface as interior zwz values have been w-masked
259            !
260            CALL nonosc_z( Kmm, pt(:,:,:,jn,Kbb), ztw, zti, p2dt )      !  monotonicity algorithm
261            !
262         CASE(  4  )                               ! 4th order COMPACT
263            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )         ! 4th order compact interpolation of T at w-point
264            DO_3D( 0, 0, 0, 0, 2, jpkm1 )
265               ztw(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
266            END_3D
267            IF( ln_linssh ) THEN
268               DO_2D( 1, 1, 1, 1 )
269                  ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)     !!gm ISF & 4th COMPACT doesn't work
270               END_2D
271            ENDIF
272            !
273         END SELECT
274         !
275         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
276            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
277               &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
278         END_3D
279         !
280         ! TEMP: These changes not necessary after trd_tra is tiled
281         IF( l_trd )  THEN       ! vertical advective trend diagnostics
282            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
283               ztrdz(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
284                  &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   &
285                  &                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
286            END_3D
287
288            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
289               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
290
291               ! TODO: TO BE TILED- trd_tra
292               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz )
293
294               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = itile )     ! Revert to tile domain
295            ENDIF
296         ENDIF
297         !
298      END DO
299      !
300   END SUBROUTINE tra_adv_ubs
301
302
303   SUBROUTINE nonosc_z( Kmm, pbef, pcc, paft, p2dt )
304      !!---------------------------------------------------------------------
305      !!                    ***  ROUTINE nonosc_z  ***
306      !!     
307      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
308      !!       scheme and the before field by a nonoscillatory algorithm
309      !!
310      !! **  Method  :   ... ???
311      !!       warning : pbef and paft must be masked, but the boundaries
312      !!       conditions on the fluxes are not necessary zalezak (1979)
313      !!       drange (1995) multi-dimensional forward-in-time and upstream-
314      !!       in-space based differencing for fluid
315      !!----------------------------------------------------------------------
316      INTEGER , INTENT(in   )                         ::   Kmm    ! time level index
317      REAL(wp), INTENT(in   )                         ::   p2dt   ! tracer time-step
318      REAL(wp),                DIMENSION(jpi,jpj,jpk) ::   pbef   ! before field
319      REAL(wp), INTENT(inout), DIMENSION(ST_2D(nn_hls)    ,jpk) ::   paft   ! after field
320      REAL(wp), INTENT(inout), DIMENSION(ST_2D(nn_hls)    ,jpk) ::   pcc    ! monotonic flux in the k direction
321      !
322      INTEGER  ::   ji, jj, jk   ! dummy loop indices
323      INTEGER  ::   ikm1         ! local integer
324      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn   ! local scalars
325      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk) ::   zbetup, zbetdo         ! 3D workspace
326      !!----------------------------------------------------------------------
327      !
328      zbig  = 1.e+38_wp
329      zrtrn = 1.e-15_wp
330      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
331      !
332      ! Search local extrema
333      ! --------------------
334      !                    ! large negative value (-zbig) inside land
335      DO_3D( 0, 0, 0, 0, 1, jpk )
336         pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1.e0 - tmask(ji,jj,jk) )
337         paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1.e0 - tmask(ji,jj,jk) )
338      END_3D
339      !
340      DO jk = 1, jpkm1     ! search maximum in neighbourhood
341         ikm1 = MAX(jk-1,1)
342         DO_2D( 0, 0, 0, 0 )
343            zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
344               &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
345               &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
346         END_2D
347      END DO
348      !                    ! large positive value (+zbig) inside land
349      DO_3D( 0, 0, 0, 0, 1, jpk )
350         pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1.e0 - tmask(ji,jj,jk) )
351         paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1.e0 - tmask(ji,jj,jk) )
352      END_3D
353      !
354      DO jk = 1, jpkm1     ! search minimum in neighbourhood
355         ikm1 = MAX(jk-1,1)
356         DO_2D( 0, 0, 0, 0 )
357            zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
358               &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
359               &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
360         END_2D
361      END DO
362      !                    ! restore masked values to zero
363      DO_3D( 0, 0, 0, 0, 1, jpk )
364         pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk)
365         paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk)
366      END_3D
367      !
368      ! Positive and negative part of fluxes and beta terms
369      ! ---------------------------------------------------
370      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
371         ! positive & negative part of the flux
372         zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
373         zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
374         ! up & down beta terms
375         zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
376         zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
377         zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
378      END_3D
379      !
380      ! monotonic flux in the k direction, i.e. pcc
381      ! -------------------------------------------
382      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
383         za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
384         zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
385         zc = 0.5 * ( 1.e0 + SIGN( 1.0_wp, pcc(ji,jj,jk) ) )
386         pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
387      END_3D
388      !
389   END SUBROUTINE nonosc_z
390
391   !!======================================================================
392END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.