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

Last change on this file since 11993 was 11993, checked in by cetlod, 9 months ago

trunk : undo bad commit. Oups

  • Property svn:keywords set to Id
File size: 20.3 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 "vectopt_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, pun, pvn, pwn,          &
49      &                                                ptb, ptn, pta, kjpt, 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 pta  with the now advective tracer trends
80      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
81      !!             - htr_adv, str_adv : 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   ) ::   kit000          ! first time step index
88      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
89      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
90      INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers
91      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
92      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport components
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
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. ln_diaptr )                                               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 jj = 1, jpjm1              ! First derivative (masked gradient)
127               DO ji = 1, fs_jpim1   ! vector opt.
128                  zeeu = e2_e1u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk)
129                  zeev = e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
130                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
131                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
132               END DO
133            END DO
134            DO jj = 2, jpjm1              ! Second derivative (divergence)
135               DO ji = fs_2, fs_jpim1   ! vector opt.
136                  zcoef = 1._wp / ( 6._wp * e3t_n(ji,jj,jk) )
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 DO
140            END DO
141            !                                   
142         END DO         
143         CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1. )   ;    CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
144         !   
145         DO jk = 1, jpkm1        !==  Horizontal advective fluxes  ==!     (UBS)
146            DO jj = 1, jpjm1
147               DO ji = 1, fs_jpim1   ! vector opt.
148                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! upstream transport (x2)
149                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
150                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
151                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
152                  !                                                  ! 2nd order centered advective fluxes (x2)
153                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
154                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
155                  !                                                  ! UBS advective fluxes
156                  ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
157                  ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
158               END DO
159            END DO
160         END DO         
161         !
162         zltu(:,:,:) = pta(:,:,:,jn)      ! store the initial trends before its update
163         !
164         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
165            DO jj = 2, jpjm1
166               DO ji = fs_2, fs_jpim1   ! vector opt.
167                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                        &
168                     &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
169                     &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
170               END DO
171            END DO
172            !                                             
173         END DO
174         !
175         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
176         !                                            ! and/or in trend diagnostic (l_trd=T)
177         !               
178         IF( l_trd ) THEN                  ! trend diagnostics
179             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pun, ptn(:,:,:,jn) )
180             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) )
181         END IF
182         !     
183         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
184         IF( l_ptr )  CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
185         !                                !  heati/salt transport
186         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
187         !
188         !
189         !                       !== vertical advective trend  ==!
190         !
191         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme
192         !
193         CASE(  2  )                   ! 2nd order FCT
194            !         
195            IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
196            !
197            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
198            DO jk = 2, jpkm1                 ! Interior value (w-masked)
199               DO jj = 1, jpj
200                  DO ji = 1, jpi
201                     zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
202                     zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
203                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk)
204                  END DO
205               END DO
206            END DO
207            IF( ln_linssh ) THEN             ! top ocean value (only in linear free surface as ztw has been w-masked)
208               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
209                  DO jj = 1, jpj
210                     DO ji = 1, jpi
211                        ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface
212                     END DO
213                  END DO   
214               ELSE                                ! no cavities: only at the ocean surface
215                  ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
216               ENDIF
217            ENDIF
218            !
219            DO jk = 1, jpkm1           !* trend and after field with monotonic scheme
220               DO jj = 2, jpjm1
221                  DO ji = fs_2, fs_jpim1   ! vector opt.
222                     ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
223                     pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
224                     zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
225                  END DO
226               END DO
227            END DO
228            CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
229            !
230            !                          !*  anti-diffusive flux : high order minus low order
231            DO jk = 2, jpkm1        ! Interior value  (w-masked)
232               DO jj = 1, jpj
233                  DO ji = 1, jpi
234                     ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   &
235                        &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
236                  END DO
237               END DO
238            END DO
239            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0
240            IF( ln_linssh )   ztw(:,:, 1 ) = 0._wp       ! only ocean surface as interior zwz values have been w-masked
241            !
242            CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
243            !
244         CASE(  4  )                               ! 4th order COMPACT
245            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point
246            DO jk = 2, jpkm1
247               DO jj = 2, jpjm1
248                  DO ji = fs_2, fs_jpim1
249                     ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
250                  END DO
251               END DO
252            END DO
253            IF( ln_linssh )   ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work
254            !
255         END SELECT
256         !
257         DO jk = 1, jpkm1        !  final trend with corrected fluxes
258            DO jj = 2, jpjm1 
259               DO ji = fs_2, fs_jpim1   ! vector opt.   
260                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
261               END DO
262            END DO
263         END DO
264         !
265         IF( l_trd )  THEN       ! vertical advective trend diagnostics
266            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
267               DO jj = 2, jpjm1
268                  DO ji = fs_2, fs_jpim1   ! vector opt.
269                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk)                          &
270                        &           + ptn(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   &
271                        &                              * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
272                  END DO
273               END DO
274            END DO
275            CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv )
276         ENDIF
277         !
278      END DO
279      !
280   END SUBROUTINE tra_adv_ubs
281
282
283   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
284      !!---------------------------------------------------------------------
285      !!                    ***  ROUTINE nonosc_z  ***
286      !!     
287      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
288      !!       scheme and the before field by a nonoscillatory algorithm
289      !!
290      !! **  Method  :   ... ???
291      !!       warning : pbef and paft must be masked, but the boundaries
292      !!       conditions on the fluxes are not necessary zalezak (1979)
293      !!       drange (1995) multi-dimensional forward-in-time and upstream-
294      !!       in-space based differencing for fluid
295      !!----------------------------------------------------------------------
296      REAL(wp), INTENT(in   )                          ::   p2dt   ! tracer time-step
297      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
298      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
299      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
300      !
301      INTEGER  ::   ji, jj, jk   ! dummy loop indices
302      INTEGER  ::   ikm1         ! local integer
303      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn   ! local scalars
304      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zbetup, zbetdo     ! 3D workspace
305      !!----------------------------------------------------------------------
306      !
307      zbig  = 1.e+40_wp
308      zrtrn = 1.e-15_wp
309      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
310      !
311      ! Search local extrema
312      ! --------------------
313      !                    ! large negative value (-zbig) inside land
314      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
315      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
316      !
317      DO jk = 1, jpkm1     ! search maximum in neighbourhood
318         ikm1 = MAX(jk-1,1)
319         DO jj = 2, jpjm1
320            DO ji = fs_2, fs_jpim1   ! vector opt.
321               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
322                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
323                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
324            END DO
325         END DO
326      END DO
327      !                    ! large positive value (+zbig) inside land
328      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
329      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
330      !
331      DO jk = 1, jpkm1     ! search minimum in neighbourhood
332         ikm1 = MAX(jk-1,1)
333         DO jj = 2, jpjm1
334            DO ji = fs_2, fs_jpim1   ! vector opt.
335               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
336                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
337                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
338            END DO
339         END DO
340      END DO
341      !                    ! restore masked values to zero
342      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
343      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
344      !
345      ! Positive and negative part of fluxes and beta terms
346      ! ---------------------------------------------------
347      DO jk = 1, jpkm1
348         DO jj = 2, jpjm1
349            DO ji = fs_2, fs_jpim1   ! vector opt.
350               ! positive & negative part of the flux
351               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
352               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
353               ! up & down beta terms
354               zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
355               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
356               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
357            END DO
358         END DO
359      END DO
360      !
361      ! monotonic flux in the k direction, i.e. pcc
362      ! -------------------------------------------
363      DO jk = 2, jpkm1
364         DO jj = 2, jpjm1
365            DO ji = fs_2, fs_jpim1   ! vector opt.
366               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
367               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
368               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
369               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
370            END DO
371         END DO
372      END DO
373      !
374   END SUBROUTINE nonosc_z
375
376   !!======================================================================
377END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.