source: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90 @ 7236

Last change on this file since 7236 was 7236, checked in by timgraham, 4 years ago

All changes related to diaptr (basin heat transports and transport components)

  • Property svn:keywords set to Id
File size: 20.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 trc_oce        ! share passive tracers/Ocean variables
17   USE trd_oce        ! trends: ocean variables
18   USE traadv_fct      ! acces to routine interp_4th_cpt
19   USE trdtra         ! trends manager: tracers
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 "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
40   !! $Id$
41   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn,          &
46      &                                                ptb, ptn, pta, kjpt, kn_ubs_v )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE tra_adv_ubs  ***
49      !!                 
50      !! ** Purpose :   Compute the now trend due to the advection of tracers
51      !!      and add it to the general trend of passive tracer equations.
52      !!
53      !! ** Method  :   The 3rd order Upstream Biased Scheme (UBS) is based on an
54      !!      upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
55      !!      It is only used in the horizontal direction.
56      !!      For example the i-component of the advective fluxes are given by :
57      !!                !  e2u e3u un ( mi(Tn) - zltu(i  ) )   if un(i) >= 0
58      !!          ztu = !  or
59      !!                !  e2u e3u un ( mi(Tn) - zltu(i+1) )   if un(i) < 0
60      !!      where zltu is the second derivative of the before temperature field:
61      !!          zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
62      !!        This results in a dissipatively dominant (i.e. hyper-diffusive)
63      !!      truncation error. The overall performance of the advection scheme
64      !!      is similar to that reported in (Farrow and Stevens, 1995).
65      !!        For stability reasons, the first term of the fluxes which corresponds
66      !!      to a second order centered scheme is evaluated using the now velocity
67      !!      (centered in time) while the second term which is the diffusive part
68      !!      of the scheme, is evaluated using the before velocity (forward in time).
69      !!      Note that UBS is not positive. Do not use it on passive tracers.
70      !!                On the vertical, the advection is evaluated using a FCT scheme,
71      !!      as the UBS have been found to be too diffusive.
72      !!                kn_ubs_v argument controles whether the FCT is based on
73      !!      a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact
74      !!      scheme (kn_ubs_v=4).
75      !!
76      !! ** Action : - update pta  with the now advective tracer trends
77      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
78      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
79      !!
80      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
81      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.
82      !!----------------------------------------------------------------------
83      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
84      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
85      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
86      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
87      INTEGER                              , INTENT(in   ) ::   kn_ubs_v        ! number of tracers
88      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
89      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean transport components
90      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
91      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
92      !
93      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
94      REAL(wp) ::   ztra, zbtr, zcoef                       ! local scalars
95      REAL(wp) ::   zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk   !   -      -
96      REAL(wp) ::   zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn    !   -      -
97      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztu, ztv, zltu, zltv, zti, ztw
98      !!----------------------------------------------------------------------
99      !
100      IF( nn_timing == 1 )  CALL timing_start('tra_adv_ubs')
101      !
102      CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zltu, zltv, zti, ztw )
103      !
104      IF( kt == kit000 )  THEN
105         IF(lwp) WRITE(numout,*)
106         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme on ', cdtype
107         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
108      ENDIF
109      !
110      l_trd = .FALSE.
111      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
112      !
113      ztw (:,:, 1 ) = 0._wp      ! surface & bottom value : set to zero for all tracers
114      zltu(:,:,jpk) = 0._wp   ;   zltv(:,:,jpk) = 0._wp
115      ztw (:,:,jpk) = 0._wp   ;   zti (:,:,jpk) = 0._wp
116      !
117      !                                                          ! ===========
118      DO jn = 1, kjpt                                            ! tracer loop
119         !                                                       ! ===========
120         !                                             
121         DO jk = 1, jpkm1        !==  horizontal laplacian of before tracer ==!
122            DO jj = 1, jpjm1              ! First derivative (masked gradient)
123               DO ji = 1, fs_jpim1   ! vector opt.
124                  zeeu = e2_e1u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk)
125                  zeev = e1_e2v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk)
126                  ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
127                  ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
128               END DO
129            END DO
130            DO jj = 2, jpjm1              ! Second derivative (divergence)
131               DO ji = fs_2, fs_jpim1   ! vector opt.
132                  zcoef = 1._wp / ( 6._wp * e3t_n(ji,jj,jk) )
133                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef
134                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef
135               END DO
136            END DO
137            !                                   
138         END DO         
139         CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
140         !   
141         DO jk = 1, jpkm1        !==  Horizontal advective fluxes  ==!     (UBS)
142            DO jj = 1, jpjm1
143               DO ji = 1, fs_jpim1   ! vector opt.
144                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! upstream transport (x2)
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         
157         !
158         zltu(:,:,:) = pta(:,:,:,jn)      ! store the initial trends before its update
159         !
160         DO jk = 1, jpkm1        !==  add the horizontal advective trend  ==!
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)  ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
166               END DO
167            END DO
168            !                                             
169         END DO
170         !
171         zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)    ! Horizontal advective trend used in vertical 2nd order FCT case
172         !                                            ! and/or in trend diagnostic (l_trd=T)
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 ) CALL dia_ptr_ohst_components( jn, 'adv', ztv(:,:,:) )
180         !
181         !                       !== vertical advective trend  ==!
182         !
183         SELECT CASE( kn_ubs_v )       ! select the vertical advection scheme
184         !
185         CASE(  2  )                   ! 2nd order FCT
186            !         
187            IF( l_trd )   zltv(:,:,:) = pta(:,:,:,jn)          ! store pta if trend diag.
188            !
189            !                          !*  upstream advection with initial mass fluxes & intermediate update  ==!
190            DO jk = 2, jpkm1                 ! Interior value (w-masked)
191               DO jj = 1, jpj
192                  DO ji = 1, jpi
193                     zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
194                     zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
195                     ztw(ji,jj,jk) = 0.5_wp * (  zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn)  ) * wmask(ji,jj,jk)
196                  END DO
197               END DO
198            END DO
199            IF( ln_linssh ) THEN             ! top ocean value (only in linear free surface as ztw has been w-masked)
200               IF( ln_isfcav ) THEN                ! top of the ice-shelf cavities and at the ocean surface
201                  DO jj = 1, jpj
202                     DO ji = 1, jpi
203                        ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface
204                     END DO
205                  END DO   
206               ELSE                                ! no cavities: only at the ocean surface
207                  ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
208               ENDIF
209            ENDIF
210            !
211            DO jk = 1, jpkm1           !* trend and after field with monotonic scheme
212               DO jj = 2, jpjm1
213                  DO ji = fs_2, fs_jpim1   ! vector opt.
214                     ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
215                     pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn) +  ztak 
216                     zti(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
217                  END DO
218               END DO
219            END DO
220            CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign)
221            !
222            !                          !*  anti-diffusive flux : high order minus low order
223            DO jk = 2, jpkm1        ! Interior value  (w-masked)
224               DO jj = 1, jpj
225                  DO ji = 1, jpi
226                     ztw(ji,jj,jk) = (   0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   &
227                        &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk)
228                  END DO
229               END DO
230            END DO
231            !                                            ! top ocean value: high order == upstream  ==>>  zwz=0
232            IF( ln_linssh )   ztw(:,:, 1 ) = 0._wp       ! only ocean surface as interior zwz values have been w-masked
233            !
234            CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt )      !  monotonicity algorithm
235            !
236         CASE(  4  )                               ! 4th order COMPACT
237            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )         ! 4th order compact interpolation of T at w-point
238            DO jk = 2, jpkm1
239               DO jj = 2, jpjm1
240                  DO ji = fs_2, fs_jpim1
241                     ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
242                  END DO
243               END DO
244            END DO
245            IF( ln_linssh )   ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)     !!gm ISF & 4th COMPACT doesn't work
246            !
247         END SELECT
248         !
249         DO jk = 1, jpkm1        !  final trend with corrected fluxes
250            DO jj = 2, jpjm1 
251               DO ji = fs_2, fs_jpim1   ! vector opt.   
252                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
253               END DO
254            END DO
255         END DO
256         !
257         IF( l_trd )  THEN       ! vertical advective trend diagnostics
258            DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
259               DO jj = 2, jpjm1
260                  DO ji = fs_2, fs_jpim1   ! vector opt.
261                     zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk)                          &
262                        &           + ptn(ji,jj,jk,jn) * (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  )   &
263                        &                              * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
264                  END DO
265               END DO
266            END DO
267            CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv )
268         ENDIF
269         !
270      END DO
271      !
272      CALL wrk_dealloc( jpi,jpj,jpk,   ztu, ztv, zltu, zltv, zti, ztw )
273      !
274      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_ubs')
275      !
276   END SUBROUTINE tra_adv_ubs
277
278
279   SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
280      !!---------------------------------------------------------------------
281      !!                    ***  ROUTINE nonosc_z  ***
282      !!     
283      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
284      !!       scheme and the before field by a nonoscillatory algorithm
285      !!
286      !! **  Method  :   ... ???
287      !!       warning : pbef and paft must be masked, but the boundaries
288      !!       conditions on the fluxes are not necessary zalezak (1979)
289      !!       drange (1995) multi-dimensional forward-in-time and upstream-
290      !!       in-space based differencing for fluid
291      !!----------------------------------------------------------------------
292      REAL(wp), INTENT(in   )                          ::   p2dt   ! 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   ! local scalars
300      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo
301      !!----------------------------------------------------------------------
302      !
303      IF( nn_timing == 1 )  CALL timing_start('nonosc_z')
304      !
305      CALL wrk_alloc( jpi,jpj,jpk,   zbetup, zbetdo )
306      !
307      zbig  = 1.e+40_wp
308      zrtrn = 1.e-15_wp
309      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
310      !
311      ! Search local extrema
312      ! --------------------
313      !                    ! large negative value (-zbig) inside land
314      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
315      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
316      !
317      DO jk = 1, jpkm1     ! search maximum in neighbourhood
318         ikm1 = MAX(jk-1,1)
319         DO jj = 2, jpjm1
320            DO ji = fs_2, fs_jpim1   ! vector opt.
321               zbetup(ji,jj,jk) = MAX(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
322                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
323                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
324            END DO
325         END DO
326      END DO
327      !                    ! large positive value (+zbig) inside land
328      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
329      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
330      !
331      DO jk = 1, jpkm1     ! search minimum in neighbourhood
332         ikm1 = MAX(jk-1,1)
333         DO jj = 2, jpjm1
334            DO ji = fs_2, fs_jpim1   ! vector opt.
335               zbetdo(ji,jj,jk) = MIN(  pbef(ji  ,jj  ,jk  ), paft(ji  ,jj  ,jk  ),   &
336                  &                     pbef(ji  ,jj  ,ikm1), pbef(ji  ,jj  ,jk+1),   &
337                  &                     paft(ji  ,jj  ,ikm1), paft(ji  ,jj  ,jk+1)  )
338            END DO
339         END DO
340      END DO
341      !                    ! restore masked values to zero
342      pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
343      paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
344      !
345      ! Positive and negative part of fluxes and beta terms
346      ! ---------------------------------------------------
347      DO jk = 1, jpkm1
348         DO jj = 2, jpjm1
349            DO ji = fs_2, fs_jpim1   ! vector opt.
350               ! positive & negative part of the flux
351               zpos = MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
352               zneg = MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
353               ! up & down beta terms
354               zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
355               zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
356               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
357            END DO
358         END DO
359      END DO
360      !
361      ! monotonic flux in the k direction, i.e. pcc
362      ! -------------------------------------------
363      DO jk = 2, jpkm1
364         DO jj = 2, jpjm1
365            DO ji = fs_2, fs_jpim1   ! vector opt.
366               za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
367               zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
368               zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
369               pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
370            END DO
371         END DO
372      END DO
373      !
374      CALL wrk_dealloc( jpi,jpj,jpk,   zbetup, zbetdo )
375      !
376      IF( nn_timing == 1 )  CALL timing_stop('nonosc_z')
377      !
378   END SUBROUTINE nonosc_z
379
380   !!======================================================================
381END MODULE traadv_ubs
Note: See TracBrowser for help on using the repository browser.