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/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/OCE/TRA – NEMO

source: NEMO/branches/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/OCE/TRA/traadv_ubs.F90 @ 13409

Last change on this file since 13409 was 13409, checked in by hadcv, 3 years ago

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