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

source: branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 5825

Last change on this file since 5825 was 5825, checked in by diovino, 9 years ago
  • Property svn:keywords set to Id
File size: 18.5 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 trc_oce        ! share passive tracers/Ocean variables
17   USE trd_oce        ! trends: ocean variables
18   USE trdtra         ! trends manager: tracers
19   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient
20   USE diaptr         ! poleward transport diagnostics
21   !
22   USE lib_mpp        ! I/O library
23   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
24   USE in_out_manager ! I/O manager
25   USE wrk_nemo       ! Memory Allocation
26   USE timing         ! Timing
27   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_adv_ubs   ! routine called by traadv module
33
34   LOGICAL :: l_trd  ! flag to compute trends or not
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
47      &                                       ptb, ptn, pta, kjpt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE tra_adv_ubs  ***
50      !!                 
51      !! ** Purpose :   Compute the now trend due to the advection of tracers
52      !!      and add it to the general trend of passive tracer equations.
53      !!
54      !! ** Method  :   The upstream biased scheme (UBS) is based on a 3rd order
55      !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
56      !!      It is only used in the horizontal direction.
57      !!      For example the i-component of the advective fluxes are given by :
58      !!                !  e2u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
59      !!          ztu = !  or
60      !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
61      !!      where zltu is the second derivative of the before temperature field:
62      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
63      !!      This results in a dissipatively dominant (i.e. hyper-diffusive)
64      !!      truncation error. The overall performance of the advection scheme
65      !!      is similar to that reported in (Farrow and Stevens, 1995).
66      !!      For stability reasons, the first term of the fluxes which corresponds
67      !!      to a second order centered scheme is evaluated using the now velocity
68      !!      (centered in time) while the second term which is the diffusive part
69      !!      of the scheme, is evaluated using the before velocity (forward in time).
70      !!      Note that UBS is not positive. Do not use it on passive tracers.
71      !!                On the vertical, the advection is evaluated using a TVD scheme,
72      !!      as the UBS have been found to be too diffusive.
73      !!
74      !! ** Action : - update (pta) with the now advective tracer trends
75      !!
76      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
77      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.
78      !!----------------------------------------------------------------------
79      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
80      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
81      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
82      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
83      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
84      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport components
85      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
86      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
87      !
88      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
89      REAL(wp) ::   ztra, zbtr, zcoef                      ! local scalars
90      REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
91      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
92      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztu, ztv, zltu, zltv, zti, ztw
93      !!----------------------------------------------------------------------
94      !
95      IF( nn_timing == 1 )  CALL timing_start('tra_adv_ubs')
96      !
97      CALL wrk_alloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw )
98      !
99      IF( kt == kit000 )  THEN
100         IF(lwp) WRITE(numout,*)
101         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
102         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
103      ENDIF
104      !
105      l_trd = .FALSE.
106      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
107      !
108      !                                                          ! ===========
109      DO jn = 1, kjpt                                            ! tracer loop
110         !                                                       ! ===========
111         ! 1. Bottom value : flux set to zero
112         ! ----------------------------------
113         zltu(:,:,jpk) = 0.e0       ;      zltv(:,:,jpk) = 0.e0
114         !                                             
115         DO jk = 1, jpkm1                                 ! Horizontal slab
116            !                                   
117            !  Laplacian
118            DO jj = 1, jpjm1            ! First derivative (gradient)
119               DO ji = 1, fs_jpim1   ! vector opt.
120                  zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
121                  zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
122                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
123                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
124               END DO
125            END DO
126            DO jj = 2, jpjm1            ! Second derivative (divergence)
127               DO ji = fs_2, fs_jpim1   ! vector opt.
128                  zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
129                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
130                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
131               END DO
132            END DO
133            !                                   
134         END DO                                           ! End of slab         
135         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
136
137         !   
138         !  Horizontal advective fluxes               
139         DO jk = 1, jpkm1                                 ! Horizontal slab
140            DO jj = 1, jpjm1
141               DO ji = 1, fs_jpim1   ! vector opt.
142                  ! upstream transport (x2)
143                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
144                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
145                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
146                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
147                  ! 2nd order centered advective fluxes (x2)
148                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
149                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
150                  ! UBS advective fluxes
151                  ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
152                  ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
153               END DO
154            END DO
155         END DO                                           ! End of slab         
156
157         zltu(:,:,:) = pta(:,:,:,jn)      ! store pta trends
158
159         DO jk = 1, jpkm1                 ! Horizontal advective trends
160            DO jj = 2, jpjm1
161               DO ji = fs_2, fs_jpim1   ! vector opt.
162                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn)                        &
163                     &             - (  ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk)    &
164                     &                + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
165               END DO
166            END DO
167            !                                             
168         END DO                                           !   End of slab
169
170         ! Horizontal trend used in tra_adv_ztvd subroutine
171         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)
172
173         !               
174         IF( l_trd ) THEN                  ! trend diagnostics
175             CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pun, ptn(:,:,:,jn) )
176             CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pvn, ptn(:,:,:,jn) )
177         END IF
178         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
179         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
180            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( ztv(:,:,:) )
181            IF( jn == jp_sal )  str_adv(:) = ptr_vj( ztv(:,:,:) )
182         ENDIF
183         
184         ! TVD scheme for the vertical direction 
185         ! ----------------------
186         IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
187
188         !  Bottom value : flux set to zero
189         ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
190
191         ! Surface value
192         IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero
193         ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface
194         ENDIF
195         !  upstream advection with initial mass fluxes & intermediate update
196         ! -------------------------------------------------------------------
197         ! Interior value
198         DO jk = 2, jpkm1
199            DO jj = 1, jpj
200               DO ji = 1, jpi
201                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
202                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
203                   ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  )
204               END DO
205            END DO
206         END DO 
207         ! update and guess with monotonic sheme
208         DO jk = 1, jpkm1
209            DO jj = 2, jpjm1
210               DO ji = fs_2, fs_jpim1   ! vector opt.
211                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
212                  ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
213                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
214                  zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
215               END DO
216            END DO
217         END DO
218         !
219         CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
220
221         !  antidiffusive flux : high order minus low order
222         ztw(:,:,1) = 0.e0       ! Surface value
223         DO jk = 2, jpkm1        ! Interior value
224            DO jj = 1, jpj
225               DO ji = 1, jpi
226                  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)
227               END DO
228            END DO
229         END DO
230         !
231         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
232
233         !  final trend with corrected fluxes
234         DO jk = 1, jpkm1
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                  ! k- vertical advective trends 
239                  ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
240                  ! added to the general tracer trends
241                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
242               END DO
243            END DO
244         END DO
245
246         !  Save the final vertical advective trends
247         IF( l_trd )  THEN                        ! vertical advective trend diagnostics
248            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
249               DO jj = 2, jpjm1
250                  DO ji = fs_2, fs_jpim1   ! vector opt.
251                     zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
252                     z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
253                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
254                  END DO
255               END DO
256            END DO
257            CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv )
258         ENDIF
259         !
260      END DO
261      !
262      CALL wrk_dealloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw )
263      !
264      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_ubs')
265      !
266   END SUBROUTINE tra_adv_ubs
267
268
269   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
270      !!---------------------------------------------------------------------
271      !!                    ***  ROUTINE nonosc_z  ***
272      !!     
273      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
274      !!       scheme and the before field by a nonoscillatory algorithm
275      !!
276      !! **  Method  :   ... ???
277      !!       warning : pbef and paft must be masked, but the boundaries
278      !!       conditions on the fluxes are not necessary zalezak (1979)
279      !!       drange (1995) multi-dimensional forward-in-time and upstream-
280      !!       in-space based differencing for fluid
281      !!----------------------------------------------------------------------
282      REAL(wp), INTENT(in   )                          ::   p2dt   ! tracer time-step
283      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
284      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
285      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
286      !
287      INTEGER  ::   ji, jj, jk   ! dummy loop indices
288      INTEGER  ::   ikm1         ! local integer
289      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn   ! local scalars
290      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo
291      !!----------------------------------------------------------------------
292      !
293      IF( nn_timing == 1 )  CALL timing_start('nonosc_z')
294      !
295      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo )
296      !
297      zbig  = 1.e+40_wp
298      zrtrn = 1.e-15_wp
299      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
300
301      ! Search local extrema
302      ! --------------------
303      ! large negative value (-zbig) inside land
304      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
305      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
306      ! search maximum in neighbourhood
307      DO jk = 1, jpkm1
308         ikm1 = MAX(jk-1,1)
309         DO jj = 2, jpjm1
310            DO ji = fs_2, fs_jpim1   ! vector opt.
311               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
312                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
313                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
314            END DO
315         END DO
316      END DO
317      ! large positive value (+zbig) inside land
318      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
319      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
320      ! search minimum in neighbourhood
321      DO jk = 1, jpkm1
322         ikm1 = MAX(jk-1,1)
323         DO jj = 2, jpjm1
324            DO ji = fs_2, fs_jpim1   ! vector opt.
325               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
326                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
327                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
328            END DO
329         END DO
330      END DO
331
332      ! restore masked values to zero
333      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
334      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
335
336
337      ! 2. Positive and negative part of fluxes and beta terms
338      ! ------------------------------------------------------
339
340      DO jk = 1, jpkm1
341         DO jj = 2, jpjm1
342            DO ji = fs_2, fs_jpim1   ! vector opt.
343               ! positive & negative part of the flux
344               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
345               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
346               ! up & down beta terms
347               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / p2dt
348               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
349               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
350            END DO
351         END DO
352      END DO
353      ! monotonic flux in the k direction, i.e. pcc
354      ! -------------------------------------------
355      DO jk = 2, jpkm1
356         DO jj = 2, jpjm1
357            DO ji = fs_2, fs_jpim1   ! vector opt.
358               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
359               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
360               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
361               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
362            END DO
363         END DO
364      END DO
365      !
366      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo )
367      !
368      IF( nn_timing == 1 )  CALL timing_stop('nonosc_z')
369      !
370   END SUBROUTINE nonosc_z
371
372   !!======================================================================
373END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.