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

source: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 2830

Last change on this file since 2830 was 2830, checked in by kpedwards, 13 years ago

Updates to average physics variables for TOP substepping.
  • Property svn:keywords set to Id
File size: 19.2 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#if defined key_top
25   USE trc, ONLY: nittrc000  !get first time step for passive tracers
26#endif
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, 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 third (UBS) is order scheme based
54      !!      on an 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      !!                !  e1u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
58      !!          zwx = !  or
59      !!                !  e1u 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, as
71      !!      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 wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
79      USE oce     , ONLY:   zwx  => ua       , zwy  => va         ! (ua,va) used as workspace
80      USE wrk_nemo, ONLY:   ztu  => wrk_3d_1 , ztv  => wrk_3d_2   ! 3D workspace
81      USE wrk_nemo, ONLY:   zltu => wrk_3d_3 , zltv => wrk_3d_4   !  -      -
82      USE wrk_nemo, ONLY:   zti  => wrk_3d_5 , ztw  => wrk_3d_6   !  -      -
83      !
84      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
85      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
86      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
87      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
88      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
89      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
90      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
91      !
92      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
93      REAL(wp) ::   ztra, zbtr, zcoef, z2dtt                       ! local scalars
94      REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
95      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
96      !!----------------------------------------------------------------------
97
98      IF( wrk_in_use(3, 1,2,3,4,5,6) )THEN
99         CALL ctl_stop('tra_adv_ubs: requested workspace arrays unavailable')   ;   RETURN
100      ENDIF
101
102#if defined key_top
103      IF( kt == nit000 .OR. (kt == nittrc000 .AND. cdtype == 'TRC'))  THEN
104#else
105      IF( kt == nit000 )  THEN
106#endif
107         IF(lwp) WRITE(numout,*)
108         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
109         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
110         !
111         l_trd = .FALSE.
112         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
113      ENDIF
114      !
115      !                                                          ! ===========
116      DO jn = 1, kjpt                                            ! tracer loop
117         !                                                       ! ===========
118         ! 1. Bottom value : flux set to zero
119         ! ----------------------------------
120         zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0
121         !                                             
122         DO jk = 1, jpkm1                                 ! Horizontal slab
123            !                                   
124            !  Laplacian
125            DO jj = 1, jpjm1            ! First derivative (gradient)
126               DO ji = 1, fs_jpim1   ! vector opt.
127                  zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
128                  zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
129                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
130                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
131               END DO
132            END DO
133            DO jj = 2, jpjm1            ! Second derivative (divergence)
134               DO ji = fs_2, fs_jpim1   ! vector opt.
135                  zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
136                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
137                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
138               END DO
139            END DO
140            !                                   
141         END DO                                           ! End of slab         
142         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
143
144         !   
145         !  Horizontal advective fluxes               
146         DO jk = 1, jpkm1                                 ! Horizontal slab
147            DO jj = 1, jpjm1
148               DO ji = 1, fs_jpim1   ! vector opt.
149                  ! upstream transport
150                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
151                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
152                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
153                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
154                  ! centered scheme
155                  zcenut = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
156                  zcenvt = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
157                  ! UBS scheme
158                  zwx(ji,jj,jk) =  zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) 
159                  zwy(ji,jj,jk) =  zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) 
160               END DO
161            END DO
162         END DO                                           ! End of slab         
163
164         zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends
165
166         ! Horizontal advective trends
167         DO jk = 1, jpkm1
168            !  Tracer flux divergence at t-point added to the general trend
169            DO jj = 2, jpjm1
170               DO ji = fs_2, fs_jpim1   ! vector opt.
171                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
172                  ! horizontal advective
173                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
174                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
175                  ! add it to the general tracer trends
176                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
177               END DO
178            END DO
179            !                                             
180         END DO                                           !   End of slab
181
182         ! Horizontal trend used in tra_adv_ztvd subroutine
183         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)
184
185         ! 3. Save the horizontal advective trends for diagnostic
186         ! ------------------------------------------------------
187         !                                 ! trend diagnostics (contribution of upstream fluxes)
188         IF( l_trd ) THEN
189             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
190             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
191         END IF
192         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
193         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
194            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
195            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
196         ENDIF
197         
198         ! TVD scheme for the vertical direction 
199         ! ----------------------
200         IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
201
202         !  Bottom value : flux set to zero
203         ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
204
205         ! Surface value
206         IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero
207         ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface
208         ENDIF
209         !  upstream advection with initial mass fluxes & intermediate update
210         ! -------------------------------------------------------------------
211         ! Interior value
212         DO jk = 2, jpkm1
213            DO jj = 1, jpj
214               DO ji = 1, jpi
215                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
216                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
217                   ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  )
218               END DO
219            END DO
220         END DO 
221         ! update and guess with monotonic sheme
222         DO jk = 1, jpkm1
223            z2dtt = p2dt(jk)
224            DO jj = 2, jpjm1
225               DO ji = fs_2, fs_jpim1   ! vector opt.
226                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
227                  ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
228                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
229                  zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
230               END DO
231            END DO
232         END DO
233         !
234         CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
235
236         !  antidiffusive flux : high order minus low order
237         ztw(:,:,1) = 0.e0       ! Surface value
238         DO jk = 2, jpkm1        ! Interior value
239            DO jj = 1, jpj
240               DO ji = 1, jpi
241                  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)
242               END DO
243            END DO
244         END DO
245         !
246         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
247
248         !  final trend with corrected fluxes
249         DO jk = 1, jpkm1
250            DO jj = 2, jpjm1 
251               DO ji = fs_2, fs_jpim1   ! vector opt.   
252                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
253                  ! k- vertical advective trends 
254                  ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
255                  ! added to the general tracer trends
256                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
257               END DO
258            END DO
259         END DO
260
261         !  Save the final vertical advective trends
262         IF( l_trd )  THEN                        ! vertical advective trend diagnostics
263            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
264               DO jj = 2, jpjm1
265                  DO ji = fs_2, fs_jpim1   ! vector opt.
266                     zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
267                     z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
268                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
269                  END DO
270               END DO
271            END DO
272            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zltv )
273         ENDIF
274         !
275      ENDDO
276      !
277      IF( wrk_not_released(3, 1,2,3,4,5,6) )   CALL ctl_stop('tra_adv_ubs: failed to release workspace arrays')
278      !
279   END SUBROUTINE tra_adv_ubs
280
281
282   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
283      !!---------------------------------------------------------------------
284      !!                    ***  ROUTINE nonosc_z  ***
285      !!     
286      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
287      !!       scheme and the before field by a nonoscillatory algorithm
288      !!
289      !! **  Method  :   ... ???
290      !!       warning : pbef and paft must be masked, but the boundaries
291      !!       conditions on the fluxes are not necessary zalezak (1979)
292      !!       drange (1995) multi-dimensional forward-in-time and upstream-
293      !!       in-space based differencing for fluid
294      !!----------------------------------------------------------------------
295      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
296      USE wrk_nemo, ONLY:   zbetup => wrk_3d_1, zbetdo => wrk_3d_2   ! 3D workspace
297      !
298      REAL(wp), INTENT(in   ), DIMENSION(jpk)          ::   p2dt   ! vertical profile of tracer time-step
299      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
300      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
301      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
302      !
303      INTEGER  ::   ji, jj, jk   ! dummy loop indices
304      INTEGER  ::   ikm1         ! local integer
305      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
306      !!----------------------------------------------------------------------
307
308      IF( wrk_in_use(3, 1,2) ) THEN
309         CALL ctl_stop('nonosc_z: requested workspace arrays unavailable')   ;   RETURN
310      ENDIF
311
312      zbig  = 1.e+40_wp
313      zrtrn = 1.e-15_wp
314      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
315
316      ! Search local extrema
317      ! --------------------
318      ! large negative value (-zbig) inside land
319      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
320      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
321      ! search maximum in neighbourhood
322      DO jk = 1, jpkm1
323         ikm1 = MAX(jk-1,1)
324         DO jj = 2, jpjm1
325            DO ji = fs_2, fs_jpim1   ! vector opt.
326               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
327                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
328                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
329            END DO
330         END DO
331      END DO
332      ! large positive value (+zbig) inside land
333      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
334      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
335      ! search minimum in neighbourhood
336      DO jk = 1, jpkm1
337         ikm1 = MAX(jk-1,1)
338         DO jj = 2, jpjm1
339            DO ji = fs_2, fs_jpim1   ! vector opt.
340               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
341                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
342                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
343            END DO
344         END DO
345      END DO
346
347      ! restore masked values to zero
348      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
349      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
350
351
352      ! 2. Positive and negative part of fluxes and beta terms
353      ! ------------------------------------------------------
354
355      DO jk = 1, jpkm1
356         z2dtt = p2dt(jk)
357         DO jj = 2, jpjm1
358            DO ji = fs_2, fs_jpim1   ! vector opt.
359               ! positive & negative part of the flux
360               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
361               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
362               ! up & down beta terms
363               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
364               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
365               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
366            END DO
367         END DO
368      END DO
369      ! monotonic flux in the k direction, i.e. pcc
370      ! -------------------------------------------
371      DO jk = 2, jpkm1
372         DO jj = 2, jpjm1
373            DO ji = fs_2, fs_jpim1   ! vector opt.
374               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
375               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
376               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
377               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
378            END DO
379         END DO
380      END DO
381      !
382      IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('nonosc_z: failed to release workspace arrays')
383      !
384   END SUBROUTINE nonosc_z
385
386   !!======================================================================
387END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.