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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

  • Property svn:keywords set to Id
File size: 19.0 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, 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      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 velocity 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      !!----------------------------------------------------------------------
94
95      IF( wrk_in_use(3, 1,2,3,4,5,6) )THEN
96         CALL ctl_stop('tra_adv_ubs: requested workspace arrays unavailable')   ;   RETURN
97      ENDIF
98
99      IF( kt == nit000 )  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         !
104         l_trd = .FALSE.
105         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
106      ENDIF
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
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                  ! centered scheme
148                  zcenut = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
149                  zcenvt = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
150                  ! UBS scheme
151                  zwx(ji,jj,jk) =  zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) 
152                  zwy(ji,jj,jk) =  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         ! Horizontal advective trends
160         DO jk = 1, jpkm1
161            !  Tracer flux divergence at t-point added to the general trend
162            DO jj = 2, jpjm1
163               DO ji = fs_2, fs_jpim1   ! vector opt.
164                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
165                  ! horizontal advective
166                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
167                     &             + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
168                  ! add it to the general tracer trends
169                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
170               END DO
171            END DO
172            !                                             
173         END DO                                           !   End of slab
174
175         ! Horizontal trend used in tra_adv_ztvd subroutine
176         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)
177
178         ! 3. Save the horizontal advective trends for diagnostic
179         ! ------------------------------------------------------
180         !                                 ! trend diagnostics (contribution of upstream fluxes)
181         IF( l_trd ) THEN
182             CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
183             CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
184         END IF
185         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
186         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
187            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
188            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
189         ENDIF
190         
191         ! TVD scheme for the vertical direction 
192         ! ----------------------
193         IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
194
195         !  Bottom value : flux set to zero
196         ztw(:,:,jpk) = 0.e0   ;   zti(:,:,jpk) = 0.e0
197
198         ! Surface value
199         IF( lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! variable volume : flux set to zero
200         ELSE                ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! free constant surface
201         ENDIF
202         !  upstream advection with initial mass fluxes & intermediate update
203         ! -------------------------------------------------------------------
204         ! Interior value
205         DO jk = 2, jpkm1
206            DO jj = 1, jpj
207               DO ji = 1, jpi
208                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
209                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
210                   ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  )
211               END DO
212            END DO
213         END DO 
214         ! update and guess with monotonic sheme
215         DO jk = 1, jpkm1
216            z2dtt = p2dt(jk)
217            DO jj = 2, jpjm1
218               DO ji = fs_2, fs_jpim1   ! vector opt.
219                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
220                  ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
221                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
222                  zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
223               END DO
224            END DO
225         END DO
226         !
227         CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
228
229         !  antidiffusive flux : high order minus low order
230         ztw(:,:,1) = 0.e0       ! Surface value
231         DO jk = 2, jpkm1        ! Interior value
232            DO jj = 1, jpj
233               DO ji = 1, jpi
234                  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)
235               END DO
236            END DO
237         END DO
238         !
239         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
240
241         !  final trend with corrected fluxes
242         DO jk = 1, jpkm1
243            DO jj = 2, jpjm1 
244               DO ji = fs_2, fs_jpim1   ! vector opt.   
245                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
246                  ! k- vertical advective trends 
247                  ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
248                  ! added to the general tracer trends
249                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
250               END DO
251            END DO
252         END DO
253
254         !  Save the final vertical advective trends
255         IF( l_trd )  THEN                        ! vertical advective trend diagnostics
256            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
257               DO jj = 2, jpjm1
258                  DO ji = fs_2, fs_jpim1   ! vector opt.
259                     zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
260                     z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
261                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
262                  END DO
263               END DO
264            END DO
265            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zltv )
266         ENDIF
267         !
268      ENDDO
269      !
270      IF( wrk_not_released(3, 1,2,3,4,5,6) )   CALL ctl_stop('tra_adv_ubs: failed to release workspace arrays')
271      !
272   END SUBROUTINE tra_adv_ubs
273
274
275   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
276      !!---------------------------------------------------------------------
277      !!                    ***  ROUTINE nonosc_z  ***
278      !!     
279      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
280      !!       scheme and the before field by a nonoscillatory algorithm
281      !!
282      !! **  Method  :   ... ???
283      !!       warning : pbef and paft must be masked, but the boundaries
284      !!       conditions on the fluxes are not necessary zalezak (1979)
285      !!       drange (1995) multi-dimensional forward-in-time and upstream-
286      !!       in-space based differencing for fluid
287      !!----------------------------------------------------------------------
288      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
289      USE wrk_nemo, ONLY:   zbetup => wrk_3d_1, zbetdo => wrk_3d_2   ! 3D workspace
290      !
291      REAL(wp), INTENT(in   ), DIMENSION(jpk)          ::   p2dt   ! vertical profile of tracer time-step
292      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
293      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
294      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
295      !
296      INTEGER  ::   ji, jj, jk   ! dummy loop indices
297      INTEGER  ::   ikm1         ! local integer
298      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
299      !!----------------------------------------------------------------------
300
301      IF( wrk_in_use(3, 1,2) ) THEN
302         CALL ctl_stop('nonosc_z: requested workspace arrays unavailable')   ;   RETURN
303      ENDIF
304
305      zbig  = 1.e+40_wp
306      zrtrn = 1.e-15_wp
307      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
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      IF( wrk_not_released(3, 1,2) )   CALL ctl_stop('nonosc_z: failed to release workspace arrays')
376      !
377   END SUBROUTINE nonosc_z
378
379   !!======================================================================
380END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.