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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 2662

Last change on this file since 2662 was 2633, checked in by trackstand2, 13 years ago

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