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_fct.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 52.9 KB
Line 
1MODULE traadv_fct
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_fct  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method)
5   !!==============================================================================
6   !! History :  3.7  !  2015-09  (L. Debreu, G. Madec)  original code (inspired from traadv_tvd.F90)
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!  tra_adv_fct    : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme
11   !!  tra_adv_fct_zts: update the tracer trend with a 3D advective trends using a 2nd order FCT scheme
12   !!                   with sub-time-stepping in the vertical direction
13   !!  nonosc         : compute monotonic tracer fluxes by a non-oscillatory algorithm
14   !!  interp_4th_cpt : 4th order compact scheme for the vertical component of the advection
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and active tracers
17   USE dom_oce        ! ocean space and time domain
18   USE trc_oce        ! share passive tracers/Ocean variables
19   USE trd_oce        ! trends: ocean variables
20   USE trdtra         ! tracers trends
21   USE diaptr         ! poleward transport diagnostics
22   USE diaar5         ! AR5 diagnostics
23   USE phycst, ONLY: rau0_rcp
24   !
25   USE in_out_manager ! I/O manager
26   USE iom
27   USE lib_mpp        ! MPP library
28   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
29   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
30   USE wrk_nemo       ! Memory Allocation
31   USE timing         ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   tra_adv_fct        ! routine called by traadv.F90
37   PUBLIC   tra_adv_fct_zts    ! routine called by traadv.F90
38   PUBLIC   interp_4th_cpt     ! routine called by traadv_cen.F90
39
40   LOGICAL  ::   l_trd   ! flag to compute trends
41   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
42   LOGICAL  ::   l_hst   ! flag to compute heat/salt transport
43   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
44
45   !                                        ! tridiag solver associated indices:
46   INTEGER, PARAMETER ::   np_NH   = 0   ! Neumann homogeneous boundary condition
47   INTEGER, PARAMETER ::   np_CEN2 = 1   ! 2nd order centered  boundary condition
48
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
53   !! $Id$
54   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn,       &
59      &                                              ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v )
60      !!----------------------------------------------------------------------
61      !!                  ***  ROUTINE tra_adv_fct  ***
62      !!
63      !! **  Purpose :   Compute the now trend due to total advection of tracers
64      !!               and add it to the general trend of tracer equations
65      !!
66      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction
67      !!               (choice through the value of kn_fct)
68      !!               - on the vertical the 4th order is a compact scheme
69      !!               - corrected flux (monotonic correction)
70      !!
71      !! ** Action : - update pta  with the now advective tracer trends
72      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
73      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
74      !!----------------------------------------------------------------------
75      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
76      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
77      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
78      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
79      INTEGER                              , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4)
80      INTEGER                              , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
81      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
82      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
83      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
84      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
85      !
86      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices 
87      REAL(wp) ::   ztra                                     ! local scalar
88      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      -
89      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      -
90      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
91      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrdz, zptry
92      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d
93      !!----------------------------------------------------------------------
94      !
95      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct')
96      !
97      CALL wrk_alloc( jpi,jpj,jpk,   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
98      !
99      IF( kt == kit000 )  THEN
100         IF(lwp) WRITE(numout,*)
101         IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype
102         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
103      ENDIF
104      !
105      l_trd = .FALSE.
106      l_hst = .FALSE.
107      l_ptr = .FALSE.
108      IF( ( cdtype == 'TRA'   .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )     l_trd = .TRUE.
109      IF(   cdtype == 'TRA'   .AND. ln_diaptr )                                              l_ptr = .TRUE. 
110      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
111         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE.
112      !
113      IF( l_trd .OR. l_hst )  THEN
114         CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
115!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
116         DO jk = 1, jpk
117            DO jj = 1, jpj
118               DO ji = 1, jpi
119                  ztrdx(ji,jj,jk) = 0._wp
120                  ztrdy(ji,jj,jk) = 0._wp
121                  ztrdz(ji,jj,jk) = 0._wp
122               END DO
123            END DO
124         END DO
125      ENDIF
126      !
127      IF( l_ptr ) THEN 
128         CALL wrk_alloc( jpi, jpj, jpk, zptry )
129!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
130         DO jk = 1, jpk
131            DO jj = 1, jpj
132               DO ji = 1, jpi
133                  zptry(ji,jj,jk) = 0._wp
134               END DO
135            END DO
136         END DO
137      ENDIF
138      !                          ! surface & bottom value : flux set to zero one for all
139!$OMP PARALLEL
140!$OMP DO schedule(static) private(jj, ji)
141      DO jj = 1, jpj
142         DO ji = 1, jpi
143            zwz(ji,jj, 1 ) = 0._wp
144            zwx(ji,jj,jpk) = 0._wp
145            zwy(ji,jj,jpk) = 0._wp
146            zwz(ji,jj,jpk) = 0._wp
147         END DO
148      END DO
149!$OMP END DO NOWAIT
150!$OMP DO schedule(static) private(jk, jj, ji)
151      DO jk = 1, jpk
152         DO jj = 1, jpj
153            DO ji = 1, jpi
154               zwi(ji,jj,jk) = 0._wp
155            END DO
156         END DO
157      END DO
158!$OMP END PARALLEL
159      !
160      DO jn = 1, kjpt            !==  loop over the tracers  ==!
161         !
162         !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
163         !                    !* upstream tracer flux in the i and j direction
164!$OMP PARALLEL
165!$OMP DO schedule(static) private(jk, jj, ji, zfp_vj, zfm_vj, zfp_ui,zfm_ui)
166         DO jk = 1, jpkm1
167            DO jj = 1, jpjm1
168               DO ji = 1, fs_jpim1   ! vector opt.
169                  ! upstream scheme
170                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
171                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
172                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
173                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
174                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
175                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
176               END DO
177            END DO
178         END DO
179!$OMP END DO NOWAIT
180         !                    !* upstream tracer flux in the k direction *!
181!$OMP DO schedule(static) private(jk, jj, ji, zfp_wk, zfm_wk)
182         DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
183            DO jj = 1, jpj
184               DO ji = 1, jpi
185                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
186                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
187                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
188               END DO
189            END DO
190         END DO
191!$OMP END PARALLEL
192         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked)
193            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface
194!$OMP PARALLEL DO schedule(static) private(jj, ji)
195               DO jj = 1, jpj
196                  DO ji = 1, jpi
197                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface
198                  END DO
199               END DO   
200            ELSE                             ! no cavities: only at the ocean surface
201!$OMP PARALLEL DO schedule(static) private(jj, ji)
202               DO jj = 1, jpj
203                  DO ji = 1, jpi
204                     zwz(ji,jj,1) = pwn(ji,jj,1) * ptb(ji,jj,1,jn)
205                  END DO
206               END DO
207            ENDIF
208         ENDIF
209         !               
210!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ztra)
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                  !                             ! total intermediate advective trends
215                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
216                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
217                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
218                  !                             ! update and guess with monotonic sheme
219                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
220                  zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
221               END DO
222            END DO
223         END DO
224         CALL lbc_lnk( zwi, 'T', 1. )  ! Lateral boundary conditions on zwi  (unchanged sign)
225         !               
226         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes)
227!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
228            DO jk = 1, jpk
229               DO jj = 1, jpj
230                  DO ji = 1, jpi
231                     ztrdx(ji,jj,jk) = zwx(ji,jj,jk)
232                     ztrdy(ji,jj,jk) = zwy(ji,jj,jk)
233                     ztrdz(ji,jj,jk) = zwz(ji,jj,jk)
234                  END DO
235               END DO
236            END DO
237         END IF
238         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
239         IF( l_ptr ) THEN
240!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
241            DO jk = 1, jpk
242               DO jj = 1, jpj
243                  DO ji = 1, jpi
244                     zptry(ji,jj,jk) = zwy(ji,jj,jk)
245                  END DO
246               END DO
247            END DO
248         END IF
249         !
250         !        !==  anti-diffusive flux : high order minus low order  ==!
251         !
252         SELECT CASE( kn_fct_h )    !* horizontal anti-diffusive fluxes
253         !
254         CASE(  2  )                   !- 2nd order centered
255!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
256            DO jk = 1, jpkm1
257               DO jj = 1, jpjm1
258                  DO ji = 1, fs_jpim1   ! vector opt.
259                     zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
260                     zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
261                  END DO
262               END DO
263            END DO
264            !
265         CASE(  4  )                   !- 4th order centered
266!$OMP PARALLEL
267!$OMP DO schedule(static) private(jj, ji)
268            DO jj = 1, jpj
269               DO ji = 1, jpi
270                  zltu(ji,jj,jpk) = 0._wp            ! Bottom value : flux set to zero
271                  zltv(ji,jj,jpk) = 0._wp
272               END DO
273            END DO
274!$OMP DO schedule(static) private(jk, jj, ji)
275            DO jk = 1, jpkm1                 ! Laplacian
276               DO jj = 1, jpjm1                    ! 1st derivative (gradient)
277                  DO ji = 1, fs_jpim1   ! vector opt.
278                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
279                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
280                  END DO
281               END DO
282               DO jj = 2, jpjm1                    ! 2nd derivative * 1/ 6
283                  DO ji = fs_2, fs_jpim1   ! vector opt.
284                     zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6
285                     zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6
286                  END DO
287               END DO
288            END DO
289!$OMP END PARALLEL
290            CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
291            !
292!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zC2t_u, zC2t_v)
293            DO jk = 1, jpkm1                 ! Horizontal advective fluxes
294               DO jj = 1, jpjm1
295                  DO ji = 1, fs_jpim1   ! vector opt.
296                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points
297                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
298                     !                                                  ! C4 minus upstream advective fluxes
299                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
300                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
301                  END DO
302               END DO
303            END DO         
304            !
305         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested
306!$OMP PARALLEL
307!$OMP DO schedule(static) private(jj, ji)
308            DO jj = 1, jpj
309               DO ji = 1, jpi
310                  ztu(ji,jj,jpk) = 0._wp             ! Bottom value : flux set to zero
311                  ztv(ji,jj,jpk) = 0._wp
312               END DO
313            END DO
314!$OMP DO schedule(static) private(jk, jj, ji)
315            DO jk = 1, jpkm1                 ! 1st derivative (gradient)
316               DO jj = 1, jpjm1
317                  DO ji = 1, fs_jpim1   ! vector opt.
318                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
319                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
320                  END DO
321               END DO
322            END DO
323!$OMP END PARALLEL
324            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
325            !
326!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zC2t_u, zC2t_v, zC4t_u, zC4t_v)
327            DO jk = 1, jpkm1                 ! Horizontal advective fluxes
328               DO jj = 2, jpjm1
329                  DO ji = 2, fs_jpim1   ! vector opt.
330                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points (x2)
331                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
332                     !                                                  ! C4 interpolation of T at u- & v-points (x2)
333                     zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) )
334                     zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) )
335                     !                                                  ! C4 minus upstream advective fluxes
336                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
337                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
338                  END DO
339               END DO
340            END DO
341            !
342         END SELECT
343         !                     
344         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values)
345         !
346         CASE(  2  )                   !- 2nd order centered
347!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
348            DO jk = 2, jpkm1   
349               DO jj = 2, jpjm1
350                  DO ji = fs_2, fs_jpim1
351                     zwz(ji,jj,jk) =  (  pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   &
352                        &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
353                  END DO
354               END DO
355            END DO
356            !
357         CASE(  4  )                   !- 4th order COMPACT
358            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )   ! zwt = COMPACT interpolation of T at w-point
359!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
360            DO jk = 2, jpkm1
361               DO jj = 2, jpjm1
362                  DO ji = fs_2, fs_jpim1
363                     zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
364                  END DO
365               END DO
366            END DO
367            !
368         END SELECT
369         IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0
370!$OMP PARALLEL DO schedule(static) private(jj, ji)
371            DO jj = 1, jpj
372               DO ji = 1, jpi
373                  zwz(ji,jj,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
374               END DO
375            END DO
376         ENDIF
377         !
378         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
379         CALL lbc_lnk( zwz, 'W',  1. )
380         !
381         !        !==  monotonicity algorithm  ==!
382         !
383         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
384         !
385         !        !==  final trend with corrected fluxes  ==!
386         !
387!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
388         DO jk = 1, jpkm1
389            DO jj = 2, jpjm1
390               DO ji = fs_2, fs_jpim1   ! vector opt. 
391                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
392                     &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
393                     &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) &
394                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
395               END DO
396            END DO
397         END DO
398         !
399         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes)
400!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
401            DO jk = 1, jpk
402               DO jj = 1, jpj
403                  DO ji = 1, jpi
404                     ztrdx(ji,jj,jk) = ztrdx(ji,jj,jk) + zwx(ji,jj,jk)  ! <<< Add to previously computed
405                     ztrdy(ji,jj,jk) = ztrdy(ji,jj,jk) + zwy(ji,jj,jk)  ! <<< Add to previously computed
406                     ztrdz(ji,jj,jk) = ztrdz(ji,jj,jk) + zwz(ji,jj,jk)  ! <<< Add to previously computed
407                  END DO
408               END DO
409            END DO
410         ENDIF
411            !
412         IF( l_trd ) THEN
413            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
414            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
415            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
416            !
417         END IF
418         !                                !  heat/salt transport
419         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
420
421         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
422         IF( l_ptr ) THEN 
423!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
424            DO jk = 1, jpk
425               DO jj = 1, jpj
426                  DO ji = 1, jpi
427                     zptry(ji,jj,jk) = zptry(ji,jj,jk) + zwy(ji,jj,jk)  ! <<< Add to previously computed
428                  END DO
429               END DO
430            END DO
431            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
432         ENDIF
433         !
434      END DO                     ! end of tracer loop
435      !
436                              CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
437      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
438      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry )
439      !
440      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct')
441      !
442   END SUBROUTINE tra_adv_fct
443
444
445   SUBROUTINE tra_adv_fct_zts( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
446      &                                                  ptb, ptn, pta, kjpt, kn_fct_zts )
447      !!----------------------------------------------------------------------
448      !!                  ***  ROUTINE tra_adv_fct_zts  ***
449      !!
450      !! **  Purpose :   Compute the now trend due to total advection of
451      !!       tracers and add it to the general trend of tracer equations
452      !!
453      !! **  Method  :   TVD ZTS scheme, i.e. 2nd order centered scheme with
454      !!       corrected flux (monotonic correction). This version use sub-
455      !!       timestepping for the vertical advection which increases stability
456      !!       when vertical metrics are small.
457      !!       note: - this advection scheme needs a leap-frog time scheme
458      !!
459      !! ** Action : - update (pta) with the now advective tracer trends
460      !!             - save the trends
461      !!----------------------------------------------------------------------
462      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
463      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
464      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
465      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
466      INTEGER                              , INTENT(in   ) ::   kn_fct_zts      ! number of number of vertical sub-timesteps
467      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
468      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
469      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
470      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
471      !
472      REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection
473      REAL(wp)                                             ::   zr_p2dt         ! reciprocal of tracer timestep
474      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices 
475      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
476      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps
477      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection
478      REAL(wp) ::   ztra            ! local scalar
479      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
480      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
481      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zwx_sav , zwy_sav
482      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav
483      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdx, ztrdy, ztrdz
484      REAL(wp), POINTER, DIMENSION(:,:,:) :: zptry
485      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrs
486      !!----------------------------------------------------------------------
487      !
488      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct_zts')
489      !
490      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
491      CALL wrk_alloc( jpi,jpj,jpk,         zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
492      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
493      !
494      IF( kt == kit000 )  THEN
495         IF(lwp) WRITE(numout,*)
496         IF(lwp) WRITE(numout,*) 'tra_adv_fct_zts : 2nd order FCT scheme with ', kn_fct_zts, ' vertical sub-timestep on ', cdtype
497         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
498      ENDIF
499      !
500      l_trd = .FALSE.
501      l_hst = .FALSE.
502      l_ptr = .FALSE.
503      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
504      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
505      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
506         &                          iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) ) l_hst = .TRUE.
507      !
508      IF( l_trd .OR. l_hst )  THEN
509         CALL wrk_alloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
510         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp
511      ENDIF
512      !
513      IF( l_ptr ) THEN 
514         CALL wrk_alloc( jpi, jpj,jpk, zptry )
515         zptry(:,:,:) = 0._wp
516      ENDIF
517      zwi(:,:,:) = 0._wp
518      z_rzts = 1._wp / REAL( kn_fct_zts, wp )
519      zr_p2dt = 1._wp / p2dt
520      !
521      ! surface & Bottom value : flux set to zero for all tracers
522      zwz(:,:, 1 ) = 0._wp
523      zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp
524      zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp
525      !
526      !                                                          ! ===========
527      DO jn = 1, kjpt                                            ! tracer loop
528         !                                                       ! ===========
529         !
530         ! Upstream advection with initial mass fluxes & intermediate update
531         DO jk = 1, jpkm1        ! upstream tracer flux in the i and j direction
532            DO jj = 1, jpjm1
533               DO ji = 1, fs_jpim1   ! vector opt.
534                  ! upstream scheme
535                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
536                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
537                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
538                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
539                  zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
540                  zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
541               END DO
542            END DO
543         END DO
544         !                       ! upstream tracer flux in the k direction
545         DO jk = 2, jpkm1              ! Interior value
546            DO jj = 1, jpj
547               DO ji = 1, jpi
548                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
549                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
550                  zwz(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)
551               END DO
552            END DO
553         END DO
554         IF( ln_linssh ) THEN          ! top value : linear free surface case only (as zwz is multiplied by wmask)
555            IF( ln_isfcav ) THEN             ! ice-shelf cavities: top value
556               DO jj = 1, jpj
557                  DO ji = 1, jpi
558                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 
559                  END DO
560               END DO   
561            ELSE                             ! no cavities, surface value
562               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
563            ENDIF
564         ENDIF
565         !
566         DO jk = 1, jpkm1         ! total advective trend
567            DO jj = 2, jpjm1
568               DO ji = fs_2, fs_jpim1   ! vector opt.
569                  !                             ! total intermediate advective trends
570                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
571                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
572                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   ) * r1_e1e2t(ji,jj)
573                  !                             ! update and guess with monotonic sheme
574                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
575                  zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
576               END DO
577            END DO
578         END DO
579         !                           
580         CALL lbc_lnk( zwi, 'T', 1. )     ! Lateral boundary conditions on zwi  (unchanged sign)
581         !               
582         IF( l_trd .OR. l_hst )  THEN                ! trend diagnostics (contribution of upstream fluxes)
583            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
584         END IF
585         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
586         IF( l_ptr )  zptry(:,:,:) = zwy(:,:,:)
587
588         ! 3. anti-diffusive flux : high order minus low order
589         ! ---------------------------------------------------
590
591         DO jk = 1, jpkm1                    !* horizontal anti-diffusive fluxes
592            !
593            DO jj = 1, jpjm1
594               DO ji = 1, fs_jpim1   ! vector opt.
595                  zwx_sav(ji,jj) = zwx(ji,jj,jk)
596                  zwy_sav(ji,jj) = zwy(ji,jj,jk)
597                  !
598                  zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
599                  zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
600               END DO
601            END DO
602            !
603            DO jj = 2, jpjm1                    ! partial horizontal divergence
604               DO ji = fs_2, fs_jpim1
605                  zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
606                     &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
607               END DO
608            END DO
609            !
610            DO jj = 1, jpjm1
611               DO ji = 1, fs_jpim1   ! vector opt.
612                  zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj)
613                  zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj)
614               END DO
615            END DO
616         END DO
617         !
618         !                                !* vertical anti-diffusive flux
619         zwz_sav(:,:,:)   = zwz(:,:,:)
620         ztrs   (:,:,:,1) = ptb(:,:,:,jn)
621         ztrs   (:,:,1,2) = ptb(:,:,1,jn)
622         ztrs   (:,:,1,3) = ptb(:,:,1,jn)
623         zwzts  (:,:,:)   = 0._wp
624         !
625         DO jl = 1, kn_fct_zts                  ! Start of sub timestepping loop
626            !
627            IF( jl == 1 ) THEN                        ! Euler forward to kick things off
628               jtb = 1   ;   jtn = 1   ;   jta = 2
629               zts(:) = p2dt * z_rzts
630               jtaken = MOD( kn_fct_zts + 1 , 2)            ! Toggle to collect every second flux
631               !                                            ! starting at jl =1 if kn_fct_zts is odd;
632               !                                            ! starting at jl =2 otherwise
633            ELSEIF( jl == 2 ) THEN                    ! First leapfrog step
634               jtb = 1   ;   jtn = 2   ;   jta = 3
635               zts(:) = 2._wp * p2dt * z_rzts
636            ELSE                                      ! Shuffle pointers for subsequent leapfrog steps
637               jtb = MOD(jtb,3) + 1
638               jtn = MOD(jtn,3) + 1
639               jta = MOD(jta,3) + 1
640            ENDIF
641            DO jk = 2, jpkm1                          ! interior value
642               DO jj = 2, jpjm1
643                  DO ji = fs_2, fs_jpim1
644                     zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) * wmask(ji,jj,jk)
645                     IF( jtaken == 0 )   zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk) * zts(jk)    ! Accumulate time-weighted vertcal flux
646                  END DO
647               END DO
648            END DO
649            IF( ln_linssh ) THEN                    ! top value (only in linear free surface case)
650               IF( ln_isfcav ) THEN                      ! ice-shelf cavities
651                  DO jj = 1, jpj
652                     DO ji = 1, jpi
653                        zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface
654                     END DO
655                  END DO   
656               ELSE                                      ! no ocean cavities
657                  zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
658               ENDIF
659            ENDIF
660            !
661            jtaken = MOD( jtaken + 1 , 2 )
662            !
663            DO jk = 2, jpkm1                             ! total advective trends
664               DO jj = 2, jpjm1
665                  DO ji = fs_2, fs_jpim1
666                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb)                                                 &
667                        &               - zts(jk) * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
668                        &                         * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
669                  END DO
670               END DO
671            END DO
672            !
673         END DO
674
675         DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping
676            DO jj = 2, jpjm1
677               DO ji = fs_2, fs_jpim1
678                  zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk)
679               END DO
680            END DO
681         END DO
682         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
683         CALL lbc_lnk( zwz, 'W',  1. )
684
685         ! 4. monotonicity algorithm
686         ! -------------------------
687         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
688
689
690         ! 5. final trend with corrected fluxes
691         ! ------------------------------------
692         DO jk = 1, jpkm1
693            DO jj = 2, jpjm1
694               DO ji = fs_2, fs_jpim1   ! vector opt. 
695                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (   zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )       &
696                     &                                    + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   )   &
697                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
698               END DO
699            END DO
700         END DO
701
702        !
703         IF( l_trd .OR. l_hst ) THEN     ! trend diagnostics (contribution of upstream fluxes)
704            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
705            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
706            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
707         ENDIF
708            !
709         IF( l_trd ) THEN
710            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
711            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
712            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
713            !
714         END IF
715         !                                             ! heat/salt transport
716         IF( l_hst )  CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
717
718         !                                            ! "Poleward" heat and salt transports (contribution of upstream fluxes)
719         IF( l_ptr ) THEN 
720            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
721            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
722         ENDIF
723         !
724      END DO
725      !
726                              CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
727                              CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
728                              CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
729      IF( l_trd .OR. l_hst )  CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
730      IF( l_ptr )             CALL wrk_dealloc( jpi, jpj, jpk, zptry )
731      !
732      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts')
733      !
734   END SUBROUTINE tra_adv_fct_zts
735
736
737   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
738      !!---------------------------------------------------------------------
739      !!                    ***  ROUTINE nonosc  ***
740      !!     
741      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
742      !!       scheme and the before field by a nonoscillatory algorithm
743      !!
744      !! **  Method  :   ... ???
745      !!       warning : pbef and paft must be masked, but the boundaries
746      !!       conditions on the fluxes are not necessary zalezak (1979)
747      !!       drange (1995) multi-dimensional forward-in-time and upstream-
748      !!       in-space based differencing for fluid
749      !!----------------------------------------------------------------------
750      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step
751      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
752      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions
753      !
754      INTEGER  ::   ji, jj, jk   ! dummy loop indices
755      INTEGER  ::   ikm1         ! local integer
756      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
757      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
758      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
759      !!----------------------------------------------------------------------
760      !
761      IF( nn_timing == 1 )  CALL timing_start('nonosc')
762      !
763      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
764      !
765      zbig  = 1.e+40_wp
766      zrtrn = 1.e-15_wp
767
768      ! Search local extrema
769      ! --------------------
770      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
771      zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ),   &
772         &        paft * tmask - zbig * ( 1._wp - tmask )  )
773      zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ),   &
774         &        paft * tmask + zbig * ( 1._wp - tmask )  )
775
776!$OMP PARALLEL
777!$OMP DO schedule(static) private(jk, jj, ji)
778      DO jk = 1, jpk
779         DO jj = 1, jpj
780            DO ji = 1, jpi
781               zbetup(ji,jj,jk) = 0._wp
782               zbetdo(ji,jj,jk) = 0._wp
783            END DO
784         END DO
785      END DO
786!$OMP DO schedule(static) private(jk, jj, ji, ikm1, zup, zdo, zpos, zneg, zbt)
787      DO jk = 1, jpkm1
788         ikm1 = MAX(jk-1,1)
789         DO jj = 2, jpjm1
790            DO ji = fs_2, fs_jpim1   ! vector opt.
791
792               ! search maximum in neighbourhood
793               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
794                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
795                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   &
796                  &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  )
797
798               ! search minimum in neighbourhood
799               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
800                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
801                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   &
802                  &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  )
803
804               ! positive part of the flux
805               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
806                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
807                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
808
809               ! negative part of the flux
810               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
811                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
812                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
813
814               ! up & down beta terms
815               zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
816               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
817               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
818            END DO
819         END DO
820      END DO
821!$OMP END PARALLEL
822      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
823
824      ! 3. monotonic flux in the i & j direction (paa & pbb)
825      ! ----------------------------------------
826!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, za, zb, zc, zav, zbv, zcv, zau, zbu, zcu)
827      DO jk = 1, jpkm1
828         DO jj = 2, jpjm1
829            DO ji = fs_2, fs_jpim1   ! vector opt.
830               zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
831               zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
832               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
833               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
834
835               zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
836               zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
837               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
838               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
839
840      ! monotonic flux in the k direction, i.e. pcc
841      ! -------------------------------------------
842               za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
843               zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
844               zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
845               pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
846            END DO
847         END DO
848      END DO
849      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
850      !
851      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
852      !
853      IF( nn_timing == 1 )  CALL timing_stop('nonosc')
854      !
855   END SUBROUTINE nonosc
856
857
858   SUBROUTINE interp_4th_cpt_org( pt_in, pt_out )
859      !!----------------------------------------------------------------------
860      !!                  ***  ROUTINE interp_4th_cpt_org  ***
861      !!
862      !! **  Purpose :   Compute the interpolation of tracer at w-point
863      !!
864      !! **  Method  :   4th order compact interpolation
865      !!----------------------------------------------------------------------
866      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! now tracer fields
867      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! now tracer field interpolated at w-pts
868      !
869      INTEGER :: ji, jj, jk   ! dummy loop integers
870      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
871      !!----------------------------------------------------------------------
872     
873      DO jk = 3, jpkm1        !==  build the three diagonal matrix  ==!
874         DO jj = 1, jpj
875            DO ji = 1, jpi
876               zwd (ji,jj,jk) = 4._wp
877               zwi (ji,jj,jk) = 1._wp
878               zws (ji,jj,jk) = 1._wp
879               zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
880               !
881               IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom
882                  zwd (ji,jj,jk) = 1._wp
883                  zwi (ji,jj,jk) = 0._wp
884                  zws (ji,jj,jk) = 0._wp
885                  zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )   
886               ENDIF
887            END DO
888         END DO
889      END DO
890      !
891      jk = 2                                          ! Switch to second order centered at top
892      DO jj = 1, jpj
893         DO ji = 1, jpi
894            zwd (ji,jj,jk) = 1._wp
895            zwi (ji,jj,jk) = 0._wp
896            zws (ji,jj,jk) = 0._wp
897            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
898         END DO
899      END DO   
900      !
901      !                       !==  tridiagonal solve  ==!
902      DO jj = 1, jpj                ! first recurrence
903         DO ji = 1, jpi
904            zwt(ji,jj,2) = zwd(ji,jj,2)
905         END DO
906      END DO
907      DO jk = 3, jpkm1
908         DO jj = 1, jpj
909            DO ji = 1, jpi
910               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
911            END DO
912         END DO
913      END DO
914      !
915      DO jj = 1, jpj                ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
916         DO ji = 1, jpi
917            pt_out(ji,jj,2) = zwrm(ji,jj,2)
918         END DO
919      END DO
920      DO jk = 3, jpkm1
921         DO jj = 1, jpj
922            DO ji = 1, jpi
923               pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
924            END DO
925         END DO
926      END DO
927
928      DO jj = 1, jpj                ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
929         DO ji = 1, jpi
930            pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
931         END DO
932      END DO
933      DO jk = jpk-2, 2, -1
934         DO jj = 1, jpj
935            DO ji = 1, jpi
936               pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
937            END DO
938         END DO
939      END DO
940      !   
941   END SUBROUTINE interp_4th_cpt_org
942   
943
944   SUBROUTINE interp_4th_cpt( pt_in, pt_out )
945      !!----------------------------------------------------------------------
946      !!                  ***  ROUTINE interp_4th_cpt  ***
947      !!
948      !! **  Purpose :   Compute the interpolation of tracer at w-point
949      !!
950      !! **  Method  :   4th order compact interpolation
951      !!----------------------------------------------------------------------
952      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! field at t-point
953      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! field interpolated at w-point
954      !
955      INTEGER ::   ji, jj, jk   ! dummy loop integers
956      INTEGER ::   ikt, ikb     ! local integers
957      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
958      !!----------------------------------------------------------------------
959      !
960      !                      !==  build the three diagonal matrix & the RHS  ==!
961      !
962      DO jk = 3, jpkm1                 ! interior (from jk=3 to jpk-1)
963         DO jj = 2, jpjm1
964            DO ji = fs_2, fs_jpim1
965               zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal
966               zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal
967               zws (ji,jj,jk) =         wmask(ji,jj,jk)                         ! upper diagonal
968               zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk)                     &   ! RHS
969                  &           *       ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) )
970            END DO
971         END DO
972      END DO
973      !
974!!gm
975!      SELECT CASE( kbc )               !* boundary condition
976!      CASE( np_NH   )   ! Neumann homogeneous at top & bottom
977!      CASE( np_CEN2 )   ! 2nd order centered  at top & bottom
978!      END SELECT
979!!gm 
980      !
981      DO jj = 2, jpjm1                 ! 2nd order centered at top & bottom
982         DO ji = fs_2, fs_jpim1
983            ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point
984            ikb = mbkt(ji,jj)                !     -   above the last wet point
985            !
986            zwd (ji,jj,ikt) = 1._wp          ! top
987            zwi (ji,jj,ikt) = 0._wp
988            zws (ji,jj,ikt) = 0._wp
989            zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
990            !
991            zwd (ji,jj,ikb) = 1._wp          ! bottom
992            zwi (ji,jj,ikb) = 0._wp
993            zws (ji,jj,ikb) = 0._wp
994            zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )           
995         END DO
996      END DO   
997      !
998      !                       !==  tridiagonal solver  ==!
999      !
1000      DO jj = 2, jpjm1              !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
1001         DO ji = fs_2, fs_jpim1
1002            zwt(ji,jj,2) = zwd(ji,jj,2)
1003         END DO
1004      END DO
1005      DO jk = 3, jpkm1
1006         DO jj = 2, jpjm1
1007            DO ji = fs_2, fs_jpim1
1008               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
1009            END DO
1010         END DO
1011      END DO
1012      !
1013      DO jj = 2, jpjm1              !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
1014         DO ji = fs_2, fs_jpim1
1015            pt_out(ji,jj,2) = zwrm(ji,jj,2)
1016         END DO
1017      END DO
1018      DO jk = 3, jpkm1
1019         DO jj = 2, jpjm1
1020            DO ji = fs_2, fs_jpim1
1021               pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
1022            END DO
1023         END DO
1024      END DO
1025
1026      DO jj = 2, jpjm1              !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
1027         DO ji = fs_2, fs_jpim1
1028            pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
1029         END DO
1030      END DO
1031      DO jk = jpk-2, 2, -1
1032         DO jj = 2, jpjm1
1033            DO ji = fs_2, fs_jpim1
1034               pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
1035            END DO
1036         END DO
1037      END DO
1038      !   
1039   END SUBROUTINE interp_4th_cpt
1040
1041
1042   SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev )
1043      !!----------------------------------------------------------------------
1044      !!                  ***  ROUTINE tridia_solver  ***
1045      !!
1046      !! **  Purpose :   solve a symmetric 3diagonal system
1047      !!
1048      !! **  Method  :   solve M.t_out = RHS(t)  where M is a tri diagonal matrix ( jpk*jpk )
1049      !!     
1050      !!             ( D_1 U_1  0   0   0  )( t_1 )   ( RHS_1 )
1051      !!             ( L_2 D_2 U_2  0   0  )( t_2 )   ( RHS_2 )
1052      !!             (  0  L_3 D_3 U_3  0  )( t_3 ) = ( RHS_3 )
1053      !!             (        ...          )( ... )   ( ...  )
1054      !!             (  0   0   0  L_k D_k )( t_k )   ( RHS_k )
1055      !!     
1056      !!        M is decomposed in the product of an upper and lower triangular matrix.
1057      !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL
1058      !!        (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
1059      !!        The solution is pta.
1060      !!        The 3d array zwt is used as a work space array.
1061      !!----------------------------------------------------------------------
1062      REAL(wp),DIMENSION(:,:,:), INTENT(in   ) ::   pD, pU, PL    ! 3-diagonal matrix
1063      REAL(wp),DIMENSION(:,:,:), INTENT(in   ) ::   pRHS          ! Right-Hand-Side
1064      REAL(wp),DIMENSION(:,:,:), INTENT(  out) ::   pt_out        !!gm field at level=F(klev)
1065      INTEGER                  , INTENT(in   ) ::   klev          ! =1 pt_out at w-level
1066      !                                                           ! =0 pt at t-level
1067      INTEGER ::   ji, jj, jk   ! dummy loop integers
1068      INTEGER ::   kstart       ! local indices
1069      REAL(wp),DIMENSION(jpi,jpj,jpk) ::   zwt   ! 3D work array
1070      !!----------------------------------------------------------------------
1071      !
1072      kstart =  1  + klev
1073      !
1074      DO jj = 2, jpjm1              !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1
1075         DO ji = fs_2, fs_jpim1
1076            zwt(ji,jj,kstart) = pD(ji,jj,kstart)
1077         END DO
1078      END DO
1079      DO jk = kstart+1, jpkm1
1080         DO jj = 2, jpjm1
1081            DO ji = fs_2, fs_jpim1
1082               zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
1083            END DO
1084         END DO
1085      END DO
1086      !
1087      DO jj = 2, jpjm1              !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
1088         DO ji = fs_2, fs_jpim1
1089            pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
1090         END DO
1091      END DO
1092      DO jk = kstart+1, jpkm1
1093         DO jj = 2, jpjm1
1094            DO ji = fs_2, fs_jpim1
1095               pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
1096            END DO
1097         END DO
1098      END DO
1099
1100      DO jj = 2, jpjm1              !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
1101         DO ji = fs_2, fs_jpim1
1102            pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
1103         END DO
1104      END DO
1105      DO jk = jpk-2, kstart, -1
1106         DO jj = 2, jpjm1
1107            DO ji = fs_2, fs_jpim1
1108               pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
1109            END DO
1110         END DO
1111      END DO
1112      !
1113   END SUBROUTINE tridia_solver
1114
1115   !!======================================================================
1116END MODULE traadv_fct
Note: See TracBrowser for help on using the repository browser.