source: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs_crs.F90 @ 10104

Last change on this file since 10104 was 10104, checked in by cbricaud, 3 years ago

clean trcadv_crs.F90 and add ubs scheme for coarsenig

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