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/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 3116

Last change on this file since 3116 was 3116, checked in by cetlod, 13 years ago

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