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 branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 3878

Last change on this file since 3878 was 3876, checked in by gm, 11 years ago

dev_r3858_CNRS3_Ediag: #927 phasing with 2011/dev_r3309_LOCEAN12_Ediag branche + mxl diag update
  • Property svn:keywords set to Id
File size: 19.1 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 trdtra         ! trends manager: tracers
19   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient
20   USE diaptr         ! poleward transport diagnostics
21   !
22   USE lib_mpp        ! I/O library
23   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
24   USE in_out_manager ! I/O manager
25   USE wrk_nemo       ! Memory Allocation
26   USE timing         ! Timing
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 or not
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
47      &                                       ptb, ptn, pta, kjpt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE tra_adv_ubs  ***
50      !!                 
51      !! ** Purpose :   Compute the now trend due to the advection of tracers
52      !!      and add it to the general trend of passive tracer equations.
53      !!
54      !! ** Method  :   The upstream biased scheme (UBS) is based on a 3rd order
55      !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
56      !!      It is only used in the horizontal direction.
57      !!      For example the i-component of the advective fluxes are given by :
58      !!                !  e2u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
59      !!          zwx = !  or
60      !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
61      !!      where zltu is the second derivative of the before temperature field:
62      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
63      !!      This results in a dissipatively dominant (i.e. hyper-diffusive)
64      !!      truncation error. The overall performance of the advection scheme
65      !!      is similar to that reported in (Farrow and Stevens, 1995).
66      !!      For stability reasons, the first term of the fluxes which corresponds
67      !!      to a second order centered scheme is evaluated using the now velocity
68      !!      (centered in time) while the second term which is the diffusive part
69      !!      of the scheme, is evaluated using the before velocity (forward in time).
70      !!      Note that UBS is not positive. Do not use it on passive tracers.
71      !!                On the vertical, the advection is evaluated using a TVD scheme,
72      !!      as the UBS have been found to be too diffusive.
73      !!
74      !! ** Action : - update (pta) with the now advective tracer trends
75      !!
76      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
77      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.
78      !!----------------------------------------------------------------------
79      USE oce     , ONLY:   zwx  => ua       , zwy  => va         ! (ua,va) used as workspace
80      !
81      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
82      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
83      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
84      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
85      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
86      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport components
87      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
88      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
89      !
90      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
91      REAL(wp) ::   ztra, zbtr, zcoef, z2dtt                       ! local scalars
92      REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
93      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
94      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztu, ztv, zltu, zltv, zti, ztw
95      !!----------------------------------------------------------------------
96      !
97      IF( nn_timing == 1 )  CALL timing_start('tra_adv_ubs')
98      !
99      CALL wrk_alloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw )
100      !
101
102      IF( kt == kit000 )  THEN
103         IF(lwp) WRITE(numout,*)
104         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
105         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
106         !
107         l_trd = .FALSE.
108         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
109      ENDIF
110      !
111      !                                                          ! ===========
112      DO jn = 1, kjpt                                            ! tracer loop
113         !                                                       ! ===========
114         ! 1. Bottom value : flux set to zero
115         ! ----------------------------------
116         zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0
117         !                                             
118         DO jk = 1, jpkm1                                 ! Horizontal slab
119            !                                   
120            !  Laplacian
121            DO jj = 1, jpjm1            ! First derivative (gradient)
122               DO ji = 1, fs_jpim1   ! vector opt.
123                  zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
124                  zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
125                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
126                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
127               END DO
128            END DO
129            DO jj = 2, jpjm1            ! Second derivative (divergence)
130               DO ji = fs_2, fs_jpim1   ! vector opt.
131                  zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
132                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
133                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
134               END DO
135            END DO
136            !                                   
137         END DO                                           ! End of slab         
138         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
139
140         !   
141         !  Horizontal advective fluxes               
142         DO jk = 1, jpkm1                                 ! Horizontal slab
143            DO jj = 1, jpjm1
144               DO ji = 1, fs_jpim1   ! vector opt.
145                  ! upstream transport (x2)
146                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
147                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
148                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
149                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
150                  ! 2nd order centered advective fluxes (x2)
151                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
152                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
153                  ! UBS advective fluxes
154                  zwx(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
155                  zwy(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
156               END DO
157            END DO
158         END DO                                           ! End of slab         
159
160         zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends
161
162         ! Horizontal advective trends
163         DO jk = 1, jpkm1
164            !  Tracer flux divergence at t-point added to the general trend
165            DO jj = 2, jpjm1
166               DO ji = fs_2, fs_jpim1   ! vector opt.
167                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
168                  ! horizontal advective
169                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
170                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
171                  ! add it to the general tracer trends
172                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
173               END DO
174            END DO
175            !                                             
176         END DO                                           !   End of slab
177
178         ! Horizontal trend used in tra_adv_ztvd subroutine
179         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)
180
181         ! 3. Save the horizontal advective trends for diagnostic
182         ! ------------------------------------------------------
183         !                                 ! trend diagnostics (contribution of upstream fluxes)
184         IF( l_trd ) THEN
185             CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
186             CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
187         END IF
188         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
189         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
190            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
191            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
192         ENDIF
193         
194         ! TVD scheme for the vertical direction 
195         ! ----------------------
196         IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
197
198         !  Bottom value : flux set to zero
199         ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
200
201         ! Surface value
202         IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero
203         ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface
204         ENDIF
205         !  upstream advection with initial mass fluxes & intermediate update
206         ! -------------------------------------------------------------------
207         ! Interior value
208         DO jk = 2, jpkm1
209            DO jj = 1, jpj
210               DO ji = 1, jpi
211                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
212                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
213                   ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  )
214               END DO
215            END DO
216         END DO 
217         ! update and guess with monotonic sheme
218         DO jk = 1, jpkm1
219            z2dtt = p2dt(jk)
220            DO jj = 2, jpjm1
221               DO ji = fs_2, fs_jpim1   ! vector opt.
222                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
223                  ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
224                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
225                  zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
226               END DO
227            END DO
228         END DO
229         !
230         CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
231
232         !  antidiffusive flux : high order minus low order
233         ztw(:,:,1) = 0.e0       ! Surface value
234         DO jk = 2, jpkm1        ! Interior value
235            DO jj = 1, jpj
236               DO ji = 1, jpi
237                  ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk)
238               END DO
239            END DO
240         END DO
241         !
242         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
243
244         !  final trend with corrected fluxes
245         DO jk = 1, jpkm1
246            DO jj = 2, jpjm1 
247               DO ji = fs_2, fs_jpim1   ! vector opt.   
248                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
249                  ! k- vertical advective trends 
250                  ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
251                  ! added to the general tracer trends
252                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
253               END DO
254            END DO
255         END DO
256
257         !  Save the final vertical advective trends
258         IF( l_trd )  THEN                        ! vertical advective trend diagnostics
259            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
260               DO jj = 2, jpjm1
261                  DO ji = fs_2, fs_jpim1   ! vector opt.
262                     zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
263                     z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
264                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
265                  END DO
266               END DO
267            END DO
268            CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv )
269         ENDIF
270         !
271      ENDDO
272      !
273      CALL wrk_dealloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw )
274      !
275      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_ubs')
276      !
277   END SUBROUTINE tra_adv_ubs
278
279
280   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
281      !!---------------------------------------------------------------------
282      !!                    ***  ROUTINE nonosc_z  ***
283      !!     
284      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
285      !!       scheme and the before field by a nonoscillatory algorithm
286      !!
287      !! **  Method  :   ... ???
288      !!       warning : pbef and paft must be masked, but the boundaries
289      !!       conditions on the fluxes are not necessary zalezak (1979)
290      !!       drange (1995) multi-dimensional forward-in-time and upstream-
291      !!       in-space based differencing for fluid
292      !!----------------------------------------------------------------------
293      !
294      REAL(wp), INTENT(in   ), DIMENSION(jpk)          ::   p2dt   ! vertical profile of tracer time-step
295      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
296      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
297      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
298      !
299      INTEGER  ::   ji, jj, jk   ! dummy loop indices
300      INTEGER  ::   ikm1         ! local integer
301      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
302      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo
303      !!----------------------------------------------------------------------
304      !
305      IF( nn_timing == 1 )  CALL timing_start('nonosc_z')
306      !
307      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo )
308      !
309
310      zbig  = 1.e+40_wp
311      zrtrn = 1.e-15_wp
312      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
313
314      ! Search local extrema
315      ! --------------------
316      ! large negative value (-zbig) inside land
317      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
318      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
319      ! search maximum in neighbourhood
320      DO jk = 1, jpkm1
321         ikm1 = MAX(jk-1,1)
322         DO jj = 2, jpjm1
323            DO ji = fs_2, fs_jpim1   ! vector opt.
324               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
325                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
326                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
327            END DO
328         END DO
329      END DO
330      ! large positive value (+zbig) inside land
331      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
332      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
333      ! search minimum in neighbourhood
334      DO jk = 1, jpkm1
335         ikm1 = MAX(jk-1,1)
336         DO jj = 2, jpjm1
337            DO ji = fs_2, fs_jpim1   ! vector opt.
338               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
339                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
340                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
341            END DO
342         END DO
343      END DO
344
345      ! restore masked values to zero
346      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
347      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
348
349
350      ! 2. Positive and negative part of fluxes and beta terms
351      ! ------------------------------------------------------
352
353      DO jk = 1, jpkm1
354         z2dtt = p2dt(jk)
355         DO jj = 2, jpjm1
356            DO ji = fs_2, fs_jpim1   ! vector opt.
357               ! positive & negative part of the flux
358               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
359               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
360               ! up & down beta terms
361               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
362               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
363               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
364            END DO
365         END DO
366      END DO
367      ! monotonic flux in the k direction, i.e. pcc
368      ! -------------------------------------------
369      DO jk = 2, jpkm1
370         DO jj = 2, jpjm1
371            DO ji = fs_2, fs_jpim1   ! vector opt.
372               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
373               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
374               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
375               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
376            END DO
377         END DO
378      END DO
379      !
380      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo )
381      !
382      IF( nn_timing == 1 )  CALL timing_stop('nonosc_z')
383      !
384   END SUBROUTINE nonosc_z
385
386   !!======================================================================
387END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.