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

Last change on this file since 13226 was 13226, checked in by orioltp, 3 months ago

Merging dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation into the trunk

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