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 trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

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