source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 3787

Last change on this file since 3787 was 3787, checked in by gm, 8 years ago

dev_MERGE_2012: #1049 & #1053 : bug corrections in traadv_eiv.F90 (AR5 diag) and in traadv_ubs.F90 (UBS scheme), resp.

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