source: NEMO/trunk/src/OCE/TRA/traadv_ubs.F90

Last change on this file was 14433, checked in by smasson, 3 weeks ago

trunk: merge dev_r14312_MPI_Interface into the trunk, #2598

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