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

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 2034

Last change on this file since 2034 was 2034, checked in by cetlod, 14 years ago

cosmetic changes
  • Property svn:eol-style set to native
  • 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
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   tra_adv_ubs   ! routine called by traadv module
28
29   LOGICAL :: l_trd  ! flag to compute trends or not
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
36   !! $Id$
37   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42   SUBROUTINE tra_adv_ubs ( kt, cdtype, 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      !!* Module used
76      USE oce         , zwx => ua   ! use ua as workspace
77      USE oce         , zwy => va   ! use va as workspace
78      !!* Arguments
79      INTEGER         , INTENT(in   )                               ::   kt              ! ocean time-step index
80      CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
81      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
82      INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
83      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
84      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
85      !!* Local declarations
86      INTEGER  ::   ji, jj, jk, jn          ! dummy loop indices
87      REAL(wp) ::   ztra, zbtr, zcoef                  ! temporary scalars
88      REAL(wp) ::   zfp_ui, zfm_ui, zcenut  !    "         "
89      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt  !    "         "    !    "         "
90      REAL(wp) ::   z2dtt, z2                   
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   !    "              "
94      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zti, ztw                !    "              "
95      !!----------------------------------------------------------------------
96
97
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme'
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      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
108      ELSE                                        ;    z2 = 2.
109      ENDIF
110      !
111      !                                                          ! ===========
112      DO jn = 1, kjpt                                            ! tracer loop
113         !                                                       ! ===========
114         ! 1. Bottom value : flux set to zero
115         ! ----------------------------------
116         zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0
117         !                                                ! ===============
118         DO jk = 1, jpkm1                                 ! Horizontal slab
119            !                                             ! ===============
120            !  Laplacian
121            ! First derivative (gradient)
122            DO jj = 1, jpjm1
123               DO ji = 1, fs_jpim1   ! vector opt.
124                  zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
125                  zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
126                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
127                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
128               END DO
129            END DO
130            ! Second derivative (divergence)
131            DO jj = 2, jpjm1
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         !                                                ! =================
141         
142         ! Lateral boundary conditions on the laplacian (zlt)   (unchanged sgn)
143         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )
144
145         !   
146         !  Horizontal advective fluxes               
147         DO jk = 1, jpkm1 
148            DO jj = 1, jpjm1
149               DO ji = 1, fs_jpim1   ! vector opt.
150                  ! upstream transport
151                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
152                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
153                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
154                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
155                  ! centered scheme
156                  zcenut = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
157                  zcenvt = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
158                  ! UBS scheme
159                  zwx(ji,jj,jk) =  zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) 
160                  zwy(ji,jj,jk) =  zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) 
161               END DO
162            END DO
163         ENDDO
164
165         zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends
166
167         ! Horizontal advective trends
168         DO jk = 1, jpkm1
169            !  Tracer flux divergence at t-point added to the general trend
170            DO jj = 2, jpjm1
171               DO ji = fs_2, fs_jpim1   ! vector opt.
172                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
173                  ! horizontal advective
174                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
175                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
176                  ! add it to the general tracer trends
177                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
178               END DO
179            END DO
180            !                                             ! ===============
181         END DO                                           !   End of slab
182         !                                                ! ===============
183
184         ! Horizontal trend used in tra_adv_ztvd subroutine
185         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)
186
187         ! 3. Save the horizontal advective trends for diagnostic
188         ! ------------------------------------------------------
189         !                                 ! trend diagnostics (contribution of upstream fluxes)
190         IF( l_trd ) THEN
191             CALL trd_tra( kt, cdtype, jn, jpt_trd_xad, zwx, pun, ptn(:,:,:,jn) )
192             CALL trd_tra( kt, cdtype, jn, jpt_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
193         END IF
194         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
195         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
196            IF( lk_zco ) THEN
197               DO jk = 1, jpkm1
198                  DO jj = 2, jpjm1
199                     DO ji = fs_2, fs_jpim1   ! vector opt.
200                       zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)                 
201                     END DO
202                  END DO
203               END DO
204            ENDIF
205            IF( jn == jp_tem )  pht_adv(:) = ptr_vj( zwy(:,:,:) )
206            IF( jn == jp_sal )  pst_adv(:) = ptr_vj( zwy(:,:,:) )
207         ENDIF
208         
209         ! TVD scheme for the vertical direction 
210         ! ----------------------
211         IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
212
213         !  Bottom value : flux set to zero
214         ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
215
216         ! Surface value
217         IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero
218         ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface
219         ENDIF
220         !  upstream advection with initial mass fluxes & intermediate update
221         ! -------------------------------------------------------------------
222         ! Interior value
223         DO jk = 2, jpkm1
224            DO jj = 1, jpj
225               DO ji = 1, jpi
226                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
227                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
228                   ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  )
229               END DO
230            END DO
231         END DO 
232         ! update and guess with monotonic sheme
233         DO jk = 1, jpkm1
234            z2dtt = z2 * rdttra(jk)
235            DO jj = 2, jpjm1
236               DO ji = fs_2, fs_jpim1   ! vector opt.
237                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
238                  ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
239                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
240                  zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
241               END DO
242            END DO
243         END DO
244         !
245         CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
246
247         !  antidiffusive flux : high order minus low order
248         ztw(:,:,1) = 0.e0       ! Surface value
249         DO jk = 2, jpkm1        ! Interior value
250            DO jj = 1, jpj
251               DO ji = 1, jpi
252                  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)
253               END DO
254            END DO
255         END DO
256         !
257         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, z2 )      !  monotonicity algorithm
258
259         !  final trend with corrected fluxes
260         DO jk = 1, jpkm1
261            DO jj = 2, jpjm1 
262               DO ji = fs_2, fs_jpim1   ! vector opt.   
263                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
264                  ! k- vertical advective trends 
265                  ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
266                  ! added to the general tracer trends
267                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
268               END DO
269            END DO
270         END DO
271
272         !  Save the final vertical advective trends
273         IF( l_trd )  THEN                        ! vertical advective trend diagnostics
274            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
275               DO jj = 2, jpjm1
276                  DO ji = fs_2, fs_jpim1   ! vector opt.
277                     zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
278                     z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
279                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
280                  END DO
281               END DO
282            END DO
283            CALL trd_tra( kt, cdtype, jn, jpt_trd_zad, zltv )
284         ENDIF
285         !
286      ENDDO
287      !
288   END SUBROUTINE tra_adv_ubs
289
290   SUBROUTINE nonosc_z( pbef, pcc, paft, prdt )
291      !!---------------------------------------------------------------------
292      !!                    ***  ROUTINE nonosc_z  ***
293      !!     
294      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
295      !!       scheme and the before field by a nonoscillatory algorithm
296      !!
297      !! **  Method  :   ... ???
298      !!       warning : pbef and paft must be masked, but the boundaries
299      !!       conditions on the fluxes are not necessary zalezak (1979)
300      !!       drange (1995) multi-dimensional forward-in-time and upstream-
301      !!       in-space based differencing for fluid
302      !!----------------------------------------------------------------------
303      REAL(wp), INTENT(in   )                          ::   prdt   ! ???
304      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
305      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
306      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
307      !!
308      INTEGER  ::   ji, jj, jk               ! dummy loop indices
309      INTEGER  ::   ikm1
310      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
311      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
312      !!----------------------------------------------------------------------
313
314      zbig = 1.e+40
315      zrtrn = 1.e-15
316      zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0
317
318      ! Search local extrema
319      ! --------------------
320      ! large negative value (-zbig) inside land
321      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
322      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
323      ! search maximum in neighbourhood
324      DO jk = 1, jpkm1
325         ikm1 = MAX(jk-1,1)
326         DO jj = 2, jpjm1
327            DO ji = fs_2, fs_jpim1   ! vector opt.
328               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
329                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
330                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
331            END DO
332         END DO
333      END DO
334      ! large positive value (+zbig) inside land
335      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
336      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
337      ! search minimum in neighbourhood
338      DO jk = 1, jpkm1
339         ikm1 = MAX(jk-1,1)
340         DO jj = 2, jpjm1
341            DO ji = fs_2, fs_jpim1   ! vector opt.
342               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
343                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
344                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
345            END DO
346         END DO
347      END DO
348
349      ! restore masked values to zero
350      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
351      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
352
353
354      ! 2. Positive and negative part of fluxes and beta terms
355      ! ------------------------------------------------------
356
357      DO jk = 1, jpkm1
358         z2dtt = prdt * rdttra(jk)
359         DO jj = 2, jpjm1
360            DO ji = fs_2, fs_jpim1   ! vector opt.
361               ! positive & negative part of the flux
362               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
363               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
364               ! up & down beta terms
365               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
366               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
367               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
368            END DO
369         END DO
370      END DO
371      ! monotonic flux in the k direction, i.e. pcc
372      ! -------------------------------------------
373      DO jk = 2, jpkm1
374         DO jj = 2, jpjm1
375            DO ji = fs_2, fs_jpim1   ! vector opt.
376
377               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
378               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
379               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
380               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
381            END DO
382         END DO
383      END DO
384      !
385   END SUBROUTINE nonosc_z
386
387   !!======================================================================
388END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.