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 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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