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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 2326

Last change on this file since 2326 was 2287, checked in by smasson, 14 years ago

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