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

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

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