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

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

Improve the merge of TRA-TRC, see ticket #717
  • Property svn:eol-style set to native
  • 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 , LOCEAN-IPSL (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence (modipsl/doc/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      !!* Module used
77      USE oce         , zwx => ua   ! use ua as workspace
78      USE oce         , zwy => va   ! use va as workspace
79      !!* Arguments
80      INTEGER         , INTENT(in   )                               ::   kt              ! ocean time-step index
81      CHARACTER(len=3), INTENT(in   )                               ::   cdtype          ! =TRA or TRC (tracer indicator)
82      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk)       ::   pun, pvn, pwn   ! 3 ocean velocity components
83      INTEGER         , INTENT(in   )                               ::   kjpt            ! number of tracers
84      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::   p2dt            ! vertical profile of tracer time-step
85      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb, ptn        ! before and now tracer fields
86      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta           ! tracer trend
87      !!* Local declarations
88      INTEGER  ::   ji, jj, jk, jn          ! dummy loop indices
89      REAL(wp) ::   ztra, zbtr, zcoef                  ! temporary scalars
90      REAL(wp) ::   zfp_ui, zfm_ui, zcenut  !    "         "
91      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt  !    "         "    !    "         "
92      REAL(wp) ::   z2dtt                   
93      REAL(wp) ::   ztak, zfp_wk, zfm_wk    !    "         "
94      REAL(wp) ::   zeeu, zeev, z_hdivn     
95      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu, ztv, zltu , zltv   !    "              "
96      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zti, ztw                !    "              "
97      !!----------------------------------------------------------------------
98
99
100      IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) )  THEN
101         IF(lwp) WRITE(numout,*)
102         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
103         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
104         !
105         l_trd = .FALSE.
106         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
107      ENDIF
108      !
109      !                                                          ! ===========
110      DO jn = 1, kjpt                                            ! tracer loop
111         !                                                       ! ===========
112         ! 1. Bottom value : flux set to zero
113         ! ----------------------------------
114         zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0
115         !                                                ! ===============
116         DO jk = 1, jpkm1                                 ! Horizontal slab
117            !                                             ! ===============
118            !  Laplacian
119            ! First derivative (gradient)
120            DO jj = 1, jpjm1
121               DO ji = 1, fs_jpim1   ! vector opt.
122                  zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
123                  zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
124                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
125                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
126               END DO
127            END DO
128            ! Second derivative (divergence)
129            DO jj = 2, jpjm1
130               DO ji = fs_2, fs_jpim1   ! vector opt.
131                  zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
132                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
133                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
134               END DO
135            END DO
136            !                                             ! =================
137         END DO                                           !    End of slab
138         !                                                ! =================
139         
140         ! Lateral boundary conditions on the laplacian (zlt)   (unchanged sgn)
141         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )
142
143         !   
144         !  Horizontal advective fluxes               
145         DO jk = 1, jpkm1 
146            DO jj = 1, jpjm1
147               DO ji = 1, fs_jpim1   ! vector opt.
148                  ! upstream transport
149                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
150                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
151                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
152                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
153                  ! centered scheme
154                  zcenut = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
155                  zcenvt = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
156                  ! UBS scheme
157                  zwx(ji,jj,jk) =  zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) 
158                  zwy(ji,jj,jk) =  zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) 
159               END DO
160            END DO
161         ENDDO
162
163         zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends
164
165         ! Horizontal advective trends
166         DO jk = 1, jpkm1
167            !  Tracer flux divergence at t-point added to the general trend
168            DO jj = 2, jpjm1
169               DO ji = fs_2, fs_jpim1   ! vector opt.
170                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
171                  ! horizontal advective
172                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
173                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
174                  ! add it to the general tracer trends
175                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
176               END DO
177            END DO
178            !                                             ! ===============
179         END DO                                           !   End of slab
180         !                                                ! ===============
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, jpt_trd_xad, zwx, pun, ptn(:,:,:,jn) )
190             CALL trd_tra( kt, cdtype, jn, jpt_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, nf_ptr ) == 0 ) ) THEN 
194            IF( lk_zco ) THEN
195               DO jk = 1, jpkm1
196                  DO jj = 2, jpjm1
197                     DO ji = fs_2, fs_jpim1   ! vector opt.
198                       zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)                 
199                     END DO
200                  END DO
201               END DO
202            ENDIF
203            IF( jn == jp_tem )  pht_adv(:) = ptr_vj( zwy(:,:,:) )
204            IF( jn == jp_sal )  pst_adv(:) = ptr_vj( zwy(:,:,:) )
205         ENDIF
206         
207         ! TVD scheme for the vertical direction 
208         ! ----------------------
209         IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
210
211         !  Bottom value : flux set to zero
212         ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
213
214         ! Surface value
215         IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero
216         ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface
217         ENDIF
218         !  upstream advection with initial mass fluxes & intermediate update
219         ! -------------------------------------------------------------------
220         ! Interior value
221         DO jk = 2, jpkm1
222            DO jj = 1, jpj
223               DO ji = 1, jpi
224                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
225                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
226                   ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  )
227               END DO
228            END DO
229         END DO 
230         ! update and guess with monotonic sheme
231         DO jk = 1, jpkm1
232            z2dtt = p2dt(jk)
233            DO jj = 2, jpjm1
234               DO ji = fs_2, fs_jpim1   ! vector opt.
235                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
236                  ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
237                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
238                  zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
239               END DO
240            END DO
241         END DO
242         !
243         CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
244
245         !  antidiffusive flux : high order minus low order
246         ztw(:,:,1) = 0.e0       ! Surface value
247         DO jk = 2, jpkm1        ! Interior value
248            DO jj = 1, jpj
249               DO ji = 1, jpi
250                  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)
251               END DO
252            END DO
253         END DO
254         !
255         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
256
257         !  final trend with corrected fluxes
258         DO jk = 1, jpkm1
259            DO jj = 2, jpjm1 
260               DO ji = fs_2, fs_jpim1   ! vector opt.   
261                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
262                  ! k- vertical advective trends 
263                  ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
264                  ! added to the general tracer trends
265                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
266               END DO
267            END DO
268         END DO
269
270         !  Save the final vertical advective trends
271         IF( l_trd )  THEN                        ! vertical advective trend diagnostics
272            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
273               DO jj = 2, jpjm1
274                  DO ji = fs_2, fs_jpim1   ! vector opt.
275                     zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
276                     z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
277                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
278                  END DO
279               END DO
280            END DO
281            CALL trd_tra( kt, cdtype, jn, jpt_trd_zad, zltv )
282         ENDIF
283         !
284      ENDDO
285      !
286   END SUBROUTINE tra_adv_ubs
287
288   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
289      !!---------------------------------------------------------------------
290      !!                    ***  ROUTINE nonosc_z  ***
291      !!     
292      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
293      !!       scheme and the before field by a nonoscillatory algorithm
294      !!
295      !! **  Method  :   ... ???
296      !!       warning : pbef and paft must be masked, but the boundaries
297      !!       conditions on the fluxes are not necessary zalezak (1979)
298      !!       drange (1995) multi-dimensional forward-in-time and upstream-
299      !!       in-space based differencing for fluid
300      !!----------------------------------------------------------------------
301      REAL(wp), INTENT(in   ), DIMENSION(jpk)          ::   p2dt            ! vertical profile of tracer time-step
302      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
303      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
304      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
305      !!
306      INTEGER  ::   ji, jj, jk               ! dummy loop indices
307      INTEGER  ::   ikm1
308      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
309      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zbetup, zbetdo
310      !!----------------------------------------------------------------------
311
312      zbig = 1.e+40
313      zrtrn = 1.e-15
314      zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0
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   END SUBROUTINE nonosc_z
383
384   !!======================================================================
385END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.