source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 5 years ago

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

  • Property svn:keywords set to Id
File size: 18.6 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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of 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, z2dtt                       ! 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 ) THEN 
180            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( ztv(:,:,:) )
181            IF( jn == jp_sal )  str_adv(:) = ptr_sj( 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            z2dtt = p2dt(jk)
210            DO jj = 2, jpjm1
211               DO ji = fs_2, fs_jpim1   ! vector opt.
212                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
213                  ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
214                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
215                  zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
216               END DO
217            END DO
218         END DO
219         !
220         CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
221
222         !  antidiffusive flux : high order minus low order
223         ztw(:,:,1) = 0.e0       ! Surface value
224         DO jk = 2, jpkm1        ! Interior value
225            DO jj = 1, jpj
226               DO ji = 1, jpi
227                  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)
228               END DO
229            END DO
230         END DO
231         !
232         CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
233
234         !  final trend with corrected fluxes
235         DO jk = 1, jpkm1
236            DO jj = 2, jpjm1 
237               DO ji = fs_2, fs_jpim1   ! vector opt.   
238                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
239                  ! k- vertical advective trends 
240                  ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
241                  ! added to the general tracer trends
242                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
243               END DO
244            END DO
245         END DO
246
247         !  Save the final vertical advective trends
248         IF( l_trd )  THEN                        ! vertical advective trend diagnostics
249            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
250               DO jj = 2, jpjm1
251                  DO ji = fs_2, fs_jpim1   ! vector opt.
252                     zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
253                     z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr
254                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
255                  END DO
256               END DO
257            END DO
258            CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv )
259         ENDIF
260         !
261      END DO
262      !
263      CALL wrk_dealloc( jpi, jpj, jpk, ztu, ztv, zltu, zltv, zti, ztw )
264      !
265      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_ubs')
266      !
267   END SUBROUTINE tra_adv_ubs
268
269
270   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
271      !!---------------------------------------------------------------------
272      !!                    ***  ROUTINE nonosc_z  ***
273      !!     
274      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
275      !!       scheme and the before field by a nonoscillatory algorithm
276      !!
277      !! **  Method  :   ... ???
278      !!       warning : pbef and paft must be masked, but the boundaries
279      !!       conditions on the fluxes are not necessary zalezak (1979)
280      !!       drange (1995) multi-dimensional forward-in-time and upstream-
281      !!       in-space based differencing for fluid
282      !!----------------------------------------------------------------------
283      REAL(wp), INTENT(in   ), DIMENSION(jpk)          ::   p2dt   ! vertical profile of tracer time-step
284      REAL(wp),                DIMENSION (jpi,jpj,jpk) ::   pbef   ! before field
285      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   paft   ! after field
286      REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) ::   pcc    ! monotonic flux in the k direction
287      !
288      INTEGER  ::   ji, jj, jk   ! dummy loop indices
289      INTEGER  ::   ikm1         ! local integer
290      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt   ! local scalars
291      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo
292      !!----------------------------------------------------------------------
293      !
294      IF( nn_timing == 1 )  CALL timing_start('nonosc_z')
295      !
296      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo )
297      !
298      zbig  = 1.e+40_wp
299      zrtrn = 1.e-15_wp
300      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
301
302      ! Search local extrema
303      ! --------------------
304      ! large negative value (-zbig) inside land
305      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
306      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
307      ! search maximum in neighbourhood
308      DO jk = 1, jpkm1
309         ikm1 = MAX(jk-1,1)
310         DO jj = 2, jpjm1
311            DO ji = fs_2, fs_jpim1   ! vector opt.
312               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
313                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
314                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
315            END DO
316         END DO
317      END DO
318      ! large positive value (+zbig) inside land
319      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
320      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
321      ! search minimum in neighbourhood
322      DO jk = 1, jpkm1
323         ikm1 = MAX(jk-1,1)
324         DO jj = 2, jpjm1
325            DO ji = fs_2, fs_jpim1   ! vector opt.
326               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
327                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
328                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
329            END DO
330         END DO
331      END DO
332
333      ! restore masked values to zero
334      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
335      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
336
337
338      ! 2. Positive and negative part of fluxes and beta terms
339      ! ------------------------------------------------------
340
341      DO jk = 1, jpkm1
342         z2dtt = p2dt(jk)
343         DO jj = 2, jpjm1
344            DO ji = fs_2, fs_jpim1   ! vector opt.
345               ! positive & negative part of the flux
346               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
347               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
348               ! up & down beta terms
349               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
350               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
351               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
352            END DO
353         END DO
354      END DO
355      ! monotonic flux in the k direction, i.e. pcc
356      ! -------------------------------------------
357      DO jk = 2, jpkm1
358         DO jj = 2, jpjm1
359            DO ji = fs_2, fs_jpim1   ! vector opt.
360               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
361               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
362               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
363               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
364            END DO
365         END DO
366      END DO
367      !
368      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo )
369      !
370      IF( nn_timing == 1 )  CALL timing_stop('nonosc_z')
371      !
372   END SUBROUTINE nonosc_z
373
374   !!======================================================================
375END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.