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/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRA/traadv_ubs.F90 @ 13228

Last change on this file since 13228 was 13228, checked in by smasson, 4 years ago

better e3: update with trunk@13227 see #2385
  • Property svn:keywords set to Id
File size: 18.4 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      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
95      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
96      !
97      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
98      REAL(wp) ::   ztra, zbtr, zcoef                       ! local scalars
99      REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
100      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
101      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu, ztv, zltu, zltv, zti, ztw   ! 3D workspace
102      !!----------------------------------------------------------------------
103      !
104      IF( kt == kit000 )  THEN
105         IF(lwp) WRITE(numout,*)
106         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
107         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
108      ENDIF
109      !
110      l_trd = .FALSE.
111      l_hst = .FALSE.
112      l_ptr = .FALSE.
113      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
114      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) )  )   l_ptr = .TRUE.
115      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
116         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
117      !
118      ztw (:,:, 1 ) = 0._wp      ! surface & bottom value : set to zero for all tracers
119      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp
120      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp
121      !
122      !                                                          ! ===========
123      DO jn = 1, kjpt                                            ! tracer loop
124         !                                                       ! ===========
125         !
126         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==!
127            DO_2D_10_10
128               zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
129               zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
130               ztu(ji,jj,jk) = zeeu * ( pt(ji+1,jj  ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
131               ztv(ji,jj,jk) = zeev * ( pt(ji  ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
132            END_2D
133            DO_2D_00_00
134               zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) )
135               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
136               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
137            END_2D
138            !
139         END DO
140         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)
141         !
142         DO_3D_10_10( 1, jpkm1 )
143            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) )      ! upstream transport (x2)
144            zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
145            zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
146            zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
147            !                                                  ! 2nd order centered advective fluxes (x2)
148            zcenut = pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm) )
149            zcenvt = pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) )
150            !                                                  ! UBS advective fluxes
151            ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
152            ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
153         END_3D
154         !
155         zltu(:,:,:) = pt(:,:,:,jn,Krhs)      ! store the initial trends before its update
156         !
157         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
158            DO_2D_00_00
159               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs)                        &
160                  &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
161                  &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) &
162                  &                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
163            END_2D
164            !
165         END DO
166         !
167         zltu(:,:,:) = pt(:,:,:,jn,Krhs) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
168         !                                            ! and/or in trend diagnostic (l_trd=T)
169         !
170         IF( l_trd ) THEN                  ! trend diagnostics
171             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztu, pU, pt(:,:,:,jn,Kmm) )
172             CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztv, pV, pt(:,:,:,jn,Kmm) )
173         END IF
174         !
175         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
176         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
177         !                                !  heati/salt transport
178         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
179         !
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(:,:,:) = pt(:,:,:,jn,Krhs)          ! store pt(:,:,:,:,Krhs) if trend diag.
188            !
189            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
190            DO_3D_11_11( 2, jpkm1 )
191               zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
192               zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
193               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)
194            END_3D
195            IF( ln_linssh ) THEN             ! top ocean value (only in linear free surface as ztw has been w-masked)
196               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
197                  DO_2D_11_11
198                     ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface
199                  END_2D
200               ELSE                                ! no cavities: only at the ocean surface
201                  ztw(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb)
202               ENDIF
203            ENDIF
204            !
205            DO_3D_00_00( 1, jpkm1 )
206               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
207                  &     * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
208               pt(ji,jj,jk,jn,Krhs) =   pt(ji,jj,jk,jn,Krhs) +  ztak
209               zti(ji,jj,jk)    = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
210            END_3D
211            CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1.0_wp )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
212            !
213            !                          !*  anti-diffusive flux : high order minus low order
214            DO_3D_11_11( 2, jpkm1 )
215               ztw(ji,jj,jk) = (   0.5_wp * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   &
216                  &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
217            END_3D
218            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0
219            IF( ln_linssh )   ztw(:,:, 1 ) = 0._wp       ! only ocean surface as interior zwz values have been w-masked
220            !
221            CALL nonosc_z( Kmm, pt(:,:,:,jn,Kbb), ztw, zti, p2dt )      !  monotonicity algorithm
222            !
223         CASE(  4  )                               ! 4th order COMPACT
224            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )         ! 4th order compact interpolation of T at w-point
225            DO_3D_00_00( 2, jpkm1 )
226               ztw(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
227            END_3D
228            IF( ln_linssh )   ztw(:,:, 1 ) = pW(:,:,1) * pt(:,:,1,jn,Kmm)     !!gm ISF & 4th COMPACT doesn't work
229            !
230         END SELECT
231         !
232         DO_3D_00_00( 1, jpkm1 )
233            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )    &
234               &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
235         END_3D
236         !
237         IF( l_trd )  THEN       ! vertical advective trend diagnostics
238            DO_3D_00_00( 1, jpkm1 )
239               zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk)                          &
240                  &           + pt(ji,jj,jk,jn,Kmm) * (  pW(ji,jj,jk) - pW(ji,jj,jk+1)  )   &
241                  &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
242            END_3D
243            CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zltv )
244         ENDIF
245         !
246      END DO
247      !
248   END SUBROUTINE tra_adv_ubs
249
250
251   SUBROUTINE nonosc_z( Kmm, pbef, pcc, paft, p2dt )
252      !!---------------------------------------------------------------------
253      !!                    ***  ROUTINE nonosc_z  ***
254      !!
255      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
256      !!       scheme and the before field by a nonoscillatory algorithm
257      !!
258      !! **  Method  :   ... ???
259      !!       warning : pbef and paft must be masked, but the boundaries
260      !!       conditions on the fluxes are not necessary zalezak (1979)
261      !!       drange (1995) multi-dimensional forward-in-time and upstream-
262      !!       in-space based differencing for fluid
263      !!----------------------------------------------------------------------
264      INTEGER , INTENT(in   )                          ::   Kmm    ! time level index
265      REAL(wp), INTENT(in   )                          ::   p2dt   ! tracer time-step
266      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
267      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
268      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
269      !
270      INTEGER  ::   ji, jj, jk   ! dummy loop indices
271      INTEGER  ::   ikm1         ! local integer
272      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn   ! local scalars
273      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zbetup, zbetdo     ! 3D workspace
274      !!----------------------------------------------------------------------
275      !
276      zbig  = 1.e+38_wp
277      zrtrn = 1.e-15_wp
278      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
279      !
280      ! Search local extrema
281      ! --------------------
282      !                    ! large negative value (-zbig) inside land
283      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
284      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
285      !
286      DO jk = 1, jpkm1     ! search maximum in neighbourhood
287         ikm1 = MAX(jk-1,1)
288         DO_2D_00_00
289            zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
290               &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
291               &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
292         END_2D
293      END DO
294      !                    ! large positive value (+zbig) inside land
295      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
296      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
297      !
298      DO jk = 1, jpkm1     ! search minimum in neighbourhood
299         ikm1 = MAX(jk-1,1)
300         DO_2D_00_00
301            zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
302               &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
303               &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
304         END_2D
305      END DO
306      !                    ! restore masked values to zero
307      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
308      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
309      !
310      ! Positive and negative part of fluxes and beta terms
311      ! ---------------------------------------------------
312      DO_3D_00_00( 1, jpkm1 )
313         ! positive & negative part of the flux
314         zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
315         zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
316         ! up & down beta terms
317         zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
318         zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
319         zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
320      END_3D
321      !
322      ! monotonic flux in the k direction, i.e. pcc
323      ! -------------------------------------------
324      DO_3D_00_00( 2, jpkm1 )
325         za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
326         zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
327         zc = 0.5 * ( 1.e0 + SIGN( 1.0_wp, pcc(ji,jj,jk) ) )
328         pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
329      END_3D
330      !
331   END SUBROUTINE nonosc_z
332
333   !!======================================================================
334END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.