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 branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 @ 7241

Last change on this file since 7241 was 7241, checked in by timgraham, 7 years ago

Added advective 2D (i/j) heat transport fields in FCT advection scheme
  • Property svn:keywords set to Id
File size: 41.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   !
23   USE in_out_manager ! I/O manager
24   USE lib_mpp        ! MPP library
25   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
26   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
27   USE wrk_nemo       ! Memory Allocation
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   tra_adv_fct        ! routine called by traadv.F90
34   PUBLIC   tra_adv_fct_zts    ! routine called by traadv.F90
35   PUBLIC   interp_4th_cpt     ! routine called by traadv_cen.F90
36
37   LOGICAL  ::   l_trd   ! flag to compute trends
38   LOGICAL  ::   l_trans   ! flag to output vertically integrated transports
39   REAL(wp) ::   r1_6 = 1._wp / 6._wp   ! =1/6
40
41   !! * Substitutions
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn,       &
51      &                                              ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v )
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE tra_adv_fct  ***
54      !!
55      !! **  Purpose :   Compute the now trend due to total advection of tracers
56      !!               and add it to the general trend of tracer equations
57      !!
58      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction
59      !!               (choice through the value of kn_fct)
60      !!               - on the vertical the 4th order is a compact scheme
61      !!               - corrected flux (monotonic correction)
62      !!
63      !! ** Action : - update pta  with the now advective tracer trends
64      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
65      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
66      !!----------------------------------------------------------------------
67      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
68      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
69      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
70      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
71      INTEGER                              , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4)
72      INTEGER                              , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4)
73      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
74      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
75      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
76      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
77      !
78      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices 
79      REAL(wp) ::   ztra                                     ! local scalar
80      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u   !   -      -
81      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v   !   -      -
82      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
83      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdx, ztrdy, ztrd, zptry
84      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d
85      !!----------------------------------------------------------------------
86      !
87      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct')
88      !
89      CALL wrk_alloc( jpi,jpj,jpk,   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
90      !
91      IF( kt == kit000 )  THEN
92         IF(lwp) WRITE(numout,*)
93         IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype
94         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
95      ENDIF
96      !
97      l_trd = .FALSE.
98      l_trans = .FALSE.
99      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
100      IF( cdtype == 'TRA' .AND. (iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") ) ) l_trans = .TRUE.
101      !
102      IF( l_trd .OR. l_trans )  THEN
103         CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
104         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp
105      ENDIF
106      !
107      IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
108         CALL wrk_alloc( jpi, jpj, jpk, zptry )
109         zptry(:,:,:) = 0._wp
110      ENDIF
111      !                          ! surface & bottom value : flux set to zero one for all
112      zwz(:,:, 1 ) = 0._wp           
113      zwx(:,:,jpk) = 0._wp   ;   zwy(:,:,jpk) = 0._wp    ;    zwz(:,:,jpk) = 0._wp
114      !
115      zwi(:,:,:) = 0._wp       
116      !
117      DO jn = 1, kjpt            !==  loop over the tracers  ==!
118         !
119         !        !==  upstream advection with initial mass fluxes & intermediate update  ==!
120         !                    !* upstream tracer flux in the i and j direction
121         DO jk = 1, jpkm1
122            DO jj = 1, jpjm1
123               DO ji = 1, fs_jpim1   ! vector opt.
124                  ! upstream scheme
125                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
126                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
127                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
128                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
129                  zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
130                  zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
131               END DO
132            END DO
133         END DO
134         !                    !* upstream tracer flux in the k direction *!
135         DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask)
136            DO jj = 1, jpj
137               DO ji = 1, jpi
138                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
139                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
140                  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)
141               END DO
142            END DO
143         END DO
144         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked)
145            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface
146               DO jj = 1, jpj
147                  DO ji = 1, jpi
148                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface
149                  END DO
150               END DO   
151            ELSE                             ! no cavities: only at the ocean surface
152               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
153            ENDIF
154         ENDIF
155         !               
156         DO jk = 1, jpkm1     !* trend and after field with monotonic scheme
157            DO jj = 2, jpjm1
158               DO ji = fs_2, fs_jpim1   ! vector opt.
159                  !                             ! total intermediate advective trends
160                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
161                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
162                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
163                  !                             ! update and guess with monotonic sheme
164                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
165                  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)
166               END DO
167            END DO
168         END DO
169         CALL lbc_lnk( zwi, 'T', 1. )  ! Lateral boundary conditions on zwi  (unchanged sign)
170         !               
171         IF( l_trd .OR. l_trans )  THEN             ! trend diagnostics (contribution of upstream fluxes)
172            ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:)
173         END IF
174         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes)
175         IF( cdtype == 'TRA' .AND. ln_diaptr )    zptry(:,:,:) = zwy(:,:,:) 
176         !
177         !        !==  anti-diffusive flux : high order minus low order  ==!
178         !
179         SELECT CASE( kn_fct_h )    !* horizontal anti-diffusive fluxes
180         !
181         CASE(  2  )                   !- 2nd order centered
182            DO jk = 1, jpkm1
183               DO jj = 1, jpjm1
184                  DO ji = 1, fs_jpim1   ! vector opt.
185                     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)
186                     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)
187                  END DO
188               END DO
189            END DO
190            !
191         CASE(  4  )                   !- 4th order centered
192            zltu(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero
193            zltv(:,:,jpk) = 0._wp
194            DO jk = 1, jpkm1                 ! Laplacian
195               DO jj = 1, jpjm1                    ! 1st derivative (gradient)
196                  DO ji = 1, fs_jpim1   ! vector opt.
197                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
198                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
199                  END DO
200               END DO
201               DO jj = 2, jpjm1                    ! 2nd derivative * 1/ 6
202                  DO ji = fs_2, fs_jpim1   ! vector opt.
203                     zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6
204                     zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6
205                  END DO
206               END DO
207            END DO
208            CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn)
209            !
210            DO jk = 1, jpkm1                 ! Horizontal advective fluxes
211               DO jj = 1, jpjm1
212                  DO ji = 1, fs_jpim1   ! vector opt.
213                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points
214                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
215                     !                                                  ! C4 minus upstream advective fluxes
216                     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)
217                     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)
218                  END DO
219               END DO
220            END DO         
221            !
222         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested
223            ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero
224            ztv(:,:,jpk) = 0._wp
225            DO jk = 1, jpkm1                 ! 1st derivative (gradient)
226               DO jj = 1, jpjm1
227                  DO ji = 1, fs_jpim1   ! vector opt.
228                     ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
229                     ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
230                  END DO
231               END DO
232            END DO
233            CALL lbc_lnk( ztu, 'U', -1. )   ;    CALL lbc_lnk( ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn)
234            !
235            DO jk = 1, jpkm1                 ! Horizontal advective fluxes
236               DO jj = 2, jpjm1
237                  DO ji = 2, fs_jpim1   ! vector opt.
238                     zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points (x2)
239                     zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn)
240                     !                                                  ! C4 interpolation of T at u- & v-points (x2)
241                     zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) )
242                     zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) )
243                     !                                                  ! C4 minus upstream advective fluxes
244                     zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
245                     zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
246                  END DO
247               END DO
248            END DO
249            !
250         END SELECT
251         !                     
252         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values)
253         !
254         CASE(  2  )                   !- 2nd order centered
255            DO jk = 2, jpkm1   
256               DO jj = 2, jpjm1
257                  DO ji = fs_2, fs_jpim1
258                     zwz(ji,jj,jk) =  (  pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   &
259                        &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk)
260                  END DO
261               END DO
262            END DO
263            !
264         CASE(  4  )                   !- 4th order COMPACT
265            CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )   ! zwt = COMPACT interpolation of T at w-point
266            DO jk = 2, jpkm1
267               DO jj = 2, jpjm1
268                  DO ji = fs_2, fs_jpim1
269                     zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
270                  END DO
271               END DO
272            END DO
273            !
274         END SELECT
275         IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0
276            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked
277         ENDIF
278         !
279         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
280         CALL lbc_lnk( zwz, 'W',  1. )
281         !
282         !        !==  monotonicity algorithm  ==!
283         !
284         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
285         !
286         !        !==  final trend with corrected fluxes  ==!
287         !
288         DO jk = 1, jpkm1
289            DO jj = 2, jpjm1
290               DO ji = fs_2, fs_jpim1   ! vector opt. 
291                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
292                     &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
293                     &                                   + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) &
294                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
295               END DO
296            END DO
297         END DO
298         !
299         IF( l_trd .OR. l_trans ) THEN     ! trend diagnostics (contribution of upstream fluxes)
300            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
301            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
302            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
303         ENDIF
304            !
305         IF( l_trd ) THEN
306            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
307            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
308            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
309            !
310         END IF
311
312         IF( l_trans .AND. jn==jp_tem ) THEN
313            CALL wrk_alloc( jpi, jpj, z2d )
314            z2d(:,:) = 0._wp 
315            DO jk = 1, jpkm1
316               DO jj = 2, jpjm1
317                  DO ji = fs_2, fs_jpim1   ! vector opt.
318                     z2d(ji,jj) = z2d(ji,jj) + ztrdx(ji,jj,jk) 
319                  END DO
320               END DO
321            END DO
322            CALL lbc_lnk( z2d, 'U', -1. )
323            CALL iom_put( "uadv_heattr", rau0_rcp * z2d )       ! heat transport in i-direction
324              !
325            z2d(:,:) = 0._wp 
326            DO jk = 1, jpkm1
327               DO jj = 2, jpjm1
328                  DO ji = fs_2, fs_jpim1   ! vector opt.
329                     z2d(ji,jj) = z2d(ji,jj) + ztrdy(ji,jj,jk) 
330                  END DO
331               END DO
332            END DO
333            CALL lbc_lnk( z2d, 'V', -1. )
334            CALL iom_put( "vadv_heattr", rau0_rcp * z2d )       ! heat transport in j-direction
335            CALL wrk_dealloc( jpi, jpj, z2d )
336         ENDIF
337         ! "Poleward" heat and salt transports (contribution of upstream fluxes)
338         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
339            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
340            CALL dia_ptr_ohst_components( jn, 'adv', zptry(:,:,:) )
341         ENDIF
342         !
343      END DO                     ! end of tracer loop
344      !
345      CALL wrk_dealloc( jpi,jpj,jpk,    zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
346      IF( l_trd .OR. l_trans )  THEN
347         CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
348      ENDIF
349      IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL wrk_dealloc( jpi, jpj, jpk, zptry )
350      !
351      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct')
352      !
353   END SUBROUTINE tra_adv_fct
354
355
356   SUBROUTINE tra_adv_fct_zts( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
357      &                                                  ptb, ptn, pta, kjpt, kn_fct_zts )
358      !!----------------------------------------------------------------------
359      !!                  ***  ROUTINE tra_adv_fct_zts  ***
360      !!
361      !! **  Purpose :   Compute the now trend due to total advection of
362      !!       tracers and add it to the general trend of tracer equations
363      !!
364      !! **  Method  :   TVD ZTS scheme, i.e. 2nd order centered scheme with
365      !!       corrected flux (monotonic correction). This version use sub-
366      !!       timestepping for the vertical advection which increases stability
367      !!       when vertical metrics are small.
368      !!       note: - this advection scheme needs a leap-frog time scheme
369      !!
370      !! ** Action : - update (pta) with the now advective tracer trends
371      !!             - save the trends
372      !!----------------------------------------------------------------------
373      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
374      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
375      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
376      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
377      INTEGER                              , INTENT(in   ) ::   kn_fct_zts      ! number of number of vertical sub-timesteps
378      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
379      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
380      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
381      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
382      !
383      REAL(wp), DIMENSION( jpk )                           ::   zts             ! length of sub-timestep for vertical advection
384      REAL(wp)                                             ::   zr_p2dt         ! reciprocal of tracer timestep
385      INTEGER  ::   ji, jj, jk, jl, jn       ! dummy loop indices 
386      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
387      INTEGER  ::   jtaken          ! toggle for collecting appropriate fluxes from sub timesteps
388      REAL(wp) ::   z_rzts          ! Fractional length of Euler forward sub-timestep for vertical advection
389      REAL(wp) ::   ztra            ! local scalar
390      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      -
391      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !   -      -
392      REAL(wp), POINTER, DIMENSION(:,:  )   ::   zwx_sav , zwy_sav
393      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav
394      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdx, ztrdy, ztrdz
395      REAL(wp), POINTER, DIMENSION(:,:,:) :: zptry
396      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrs
397      !!----------------------------------------------------------------------
398      !
399      IF( nn_timing == 1 )  CALL timing_start('tra_adv_fct_zts')
400      !
401      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
402      CALL wrk_alloc( jpi,jpj,jpk,         zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
403      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
404      !
405      IF( kt == kit000 )  THEN
406         IF(lwp) WRITE(numout,*)
407         IF(lwp) WRITE(numout,*) 'tra_adv_fct_zts : 2nd order FCT scheme with ', kn_fct_zts, ' vertical sub-timestep on ', cdtype
408         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
409      ENDIF
410      !
411      l_trd = .FALSE.
412      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
413      !
414      IF( l_trd )  THEN
415         CALL wrk_alloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
416         ztrdx(:,:,:) = 0._wp  ;    ztrdy(:,:,:) = 0._wp  ;   ztrdz(:,:,:) = 0._wp
417      ENDIF
418      !
419      IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
420         CALL wrk_alloc( jpi, jpj,jpk, zptry )
421         zptry(:,:,:) = 0._wp
422      ENDIF
423      zwi(:,:,:) = 0._wp
424      z_rzts = 1._wp / REAL( kn_fct_zts, wp )
425      zr_p2dt = 1._wp / p2dt
426      !
427      ! surface & Bottom value : flux set to zero for all tracers
428      zwz(:,:, 1 ) = 0._wp
429      zwx(:,:,jpk) = 0._wp   ;    zwz(:,:,jpk) = 0._wp
430      zwy(:,:,jpk) = 0._wp   ;    zwi(:,:,jpk) = 0._wp
431      !
432      !                                                          ! ===========
433      DO jn = 1, kjpt                                            ! tracer loop
434         !                                                       ! ===========
435         !
436         ! Upstream advection with initial mass fluxes & intermediate update
437         DO jk = 1, jpkm1        ! upstream tracer flux in the i and j direction
438            DO jj = 1, jpjm1
439               DO ji = 1, fs_jpim1   ! vector opt.
440                  ! upstream scheme
441                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
442                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
443                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
444                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
445                  zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) )
446                  zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) )
447               END DO
448            END DO
449         END DO
450         !                       ! upstream tracer flux in the k direction
451         DO jk = 2, jpkm1              ! Interior value
452            DO jj = 1, jpj
453               DO ji = 1, jpi
454                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
455                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
456                  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)
457               END DO
458            END DO
459         END DO
460         IF( ln_linssh ) THEN          ! top value : linear free surface case only (as zwz is multiplied by wmask)
461            IF( ln_isfcav ) THEN             ! ice-shelf cavities: top value
462               DO jj = 1, jpj
463                  DO ji = 1, jpi
464                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 
465                  END DO
466               END DO   
467            ELSE                             ! no cavities, surface value
468               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
469            ENDIF
470         ENDIF
471         !
472         DO jk = 1, jpkm1         ! total advective trend
473            DO jj = 2, jpjm1
474               DO ji = fs_2, fs_jpim1   ! vector opt.
475                  !                             ! total intermediate advective trends
476                  ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
477                     &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
478                     &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   ) * r1_e1e2t(ji,jj)
479                  !                             ! update and guess with monotonic sheme
480                  pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
481                  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)
482               END DO
483            END DO
484         END DO
485         !                           
486         CALL lbc_lnk( zwi, 'T', 1. )     ! Lateral boundary conditions on zwi  (unchanged sign)
487         !               
488         IF( l_trd )  THEN                ! trend diagnostics (contribution of upstream fluxes)
489            ztrdx(:,:,:) = zwx(:,:,:)   ;    ztrdy(:,:,:) = zwy(:,:,:)  ;   ztrdz(:,:,:) = zwz(:,:,:)
490         END IF
491         !                                ! "Poleward" heat and salt transports (contribution of upstream fluxes)
492         IF( cdtype == 'TRA' .AND. ln_diaptr )  zptry(:,:,:) = zwy(:,:,:)
493
494         ! 3. anti-diffusive flux : high order minus low order
495         ! ---------------------------------------------------
496
497         DO jk = 1, jpkm1                    !* horizontal anti-diffusive fluxes
498            !
499            DO jj = 1, jpjm1
500               DO ji = 1, fs_jpim1   ! vector opt.
501                  zwx_sav(ji,jj) = zwx(ji,jj,jk)
502                  zwy_sav(ji,jj) = zwy(ji,jj,jk)
503                  !
504                  zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
505                  zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
506               END DO
507            END DO
508            !
509            DO jj = 2, jpjm1                    ! partial horizontal divergence
510               DO ji = fs_2, fs_jpim1
511                  zhdiv(ji,jj,jk) = (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
512                     &               + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
513               END DO
514            END DO
515            !
516            DO jj = 1, jpjm1
517               DO ji = 1, fs_jpim1   ! vector opt.
518                  zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj)
519                  zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj)
520               END DO
521            END DO
522         END DO
523         !
524         !                                !* vertical anti-diffusive flux
525         zwz_sav(:,:,:)   = zwz(:,:,:)
526         ztrs   (:,:,:,1) = ptb(:,:,:,jn)
527         ztrs   (:,:,1,2) = ptb(:,:,1,jn)
528         ztrs   (:,:,1,3) = ptb(:,:,1,jn)
529         zwzts  (:,:,:)   = 0._wp
530         !
531         DO jl = 1, kn_fct_zts                  ! Start of sub timestepping loop
532            !
533            IF( jl == 1 ) THEN                        ! Euler forward to kick things off
534               jtb = 1   ;   jtn = 1   ;   jta = 2
535               zts(:) = p2dt * z_rzts
536               jtaken = MOD( kn_fct_zts + 1 , 2)            ! Toggle to collect every second flux
537               !                                            ! starting at jl =1 if kn_fct_zts is odd;
538               !                                            ! starting at jl =2 otherwise
539            ELSEIF( jl == 2 ) THEN                    ! First leapfrog step
540               jtb = 1   ;   jtn = 2   ;   jta = 3
541               zts(:) = 2._wp * p2dt * z_rzts
542            ELSE                                      ! Shuffle pointers for subsequent leapfrog steps
543               jtb = MOD(jtb,3) + 1
544               jtn = MOD(jtn,3) + 1
545               jta = MOD(jta,3) + 1
546            ENDIF
547            DO jk = 2, jpkm1                          ! interior value
548               DO jj = 2, jpjm1
549                  DO ji = fs_2, fs_jpim1
550                     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)
551                     IF( jtaken == 0 )   zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk) * zts(jk)    ! Accumulate time-weighted vertcal flux
552                  END DO
553               END DO
554            END DO
555            IF( ln_linssh ) THEN                    ! top value (only in linear free surface case)
556               IF( ln_isfcav ) THEN                      ! ice-shelf cavities
557                  DO jj = 1, jpj
558                     DO ji = 1, jpi
559                        zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface
560                     END DO
561                  END DO   
562               ELSE                                      ! no ocean cavities
563                  zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
564               ENDIF
565            ENDIF
566            !
567            jtaken = MOD( jtaken + 1 , 2 )
568            !
569            DO jk = 2, jpkm1                             ! total advective trends
570               DO jj = 2, jpjm1
571                  DO ji = fs_2, fs_jpim1
572                     ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb)                                                 &
573                        &               - zts(jk) * (  zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
574                        &                         * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
575                  END DO
576               END DO
577            END DO
578            !
579         END DO
580
581         DO jk = 2, jpkm1          ! Anti-diffusive vertical flux using average flux from the sub-timestepping
582            DO jj = 2, jpjm1
583               DO ji = fs_2, fs_jpim1
584                  zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk)
585               END DO
586            END DO
587         END DO
588         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions
589         CALL lbc_lnk( zwz, 'W',  1. )
590
591         ! 4. monotonicity algorithm
592         ! -------------------------
593         CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
594
595
596         ! 5. final trend with corrected fluxes
597         ! ------------------------------------
598         DO jk = 1, jpkm1
599            DO jj = 2, jpjm1
600               DO ji = fs_2, fs_jpim1   ! vector opt. 
601                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (   zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )       &
602                     &                                    + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)   )   &
603                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
604               END DO
605            END DO
606         END DO
607
608         !                                 ! trend diagnostics (contribution of upstream fluxes)
609         IF( l_trd )  THEN
610            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< Add to previously computed
611            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  ! <<< Add to previously computed
612            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! <<< Add to previously computed
613            !
614            CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )   
615            CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
616            CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
617            !
618            CALL wrk_dealloc( jpi,jpj,jpk,   ztrdx, ztrdy, ztrdz )
619         END IF
620         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
621         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
622            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) 
623            CALL dia_ptr_ohst_components( jn, 'adv', zptry(:,:,:) )
624         ENDIF
625         !
626      END DO
627      !
628      CALL wrk_alloc( jpi,jpj,             zwx_sav, zwy_sav )
629      CALL wrk_alloc( jpi,jpj, jpk,        zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
630      CALL wrk_alloc( jpi,jpj,jpk,kjpt+1,  ztrs )
631      IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL wrk_dealloc( jpi, jpj, jpk, zptry )
632      !
633      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_fct_zts')
634      !
635   END SUBROUTINE tra_adv_fct_zts
636
637
638   SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
639      !!---------------------------------------------------------------------
640      !!                    ***  ROUTINE nonosc  ***
641      !!     
642      !! **  Purpose :   compute monotonic tracer fluxes from the upstream
643      !!       scheme and the before field by a nonoscillatory algorithm
644      !!
645      !! **  Method  :   ... ???
646      !!       warning : pbef and paft must be masked, but the boundaries
647      !!       conditions on the fluxes are not necessary zalezak (1979)
648      !!       drange (1995) multi-dimensional forward-in-time and upstream-
649      !!       in-space based differencing for fluid
650      !!----------------------------------------------------------------------
651      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step
652      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field
653      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions
654      !
655      INTEGER  ::   ji, jj, jk   ! dummy loop indices
656      INTEGER  ::   ikm1         ! local integer
657      REAL(wp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars
658      REAL(wp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      -
659      REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
660      !!----------------------------------------------------------------------
661      !
662      IF( nn_timing == 1 )  CALL timing_start('nonosc')
663      !
664      CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
665      !
666      zbig  = 1.e+40_wp
667      zrtrn = 1.e-15_wp
668      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp
669
670      ! Search local extrema
671      ! --------------------
672      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
673      zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ),   &
674         &        paft * tmask - zbig * ( 1._wp - tmask )  )
675      zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ),   &
676         &        paft * tmask + zbig * ( 1._wp - tmask )  )
677
678      DO jk = 1, jpkm1
679         ikm1 = MAX(jk-1,1)
680         DO jj = 2, jpjm1
681            DO ji = fs_2, fs_jpim1   ! vector opt.
682
683               ! search maximum in neighbourhood
684               zup = MAX(  zbup(ji  ,jj  ,jk  ),   &
685                  &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   &
686                  &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   &
687                  &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  )
688
689               ! search minimum in neighbourhood
690               zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   &
691                  &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   &
692                  &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   &
693                  &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  )
694
695               ! positive part of the flux
696               zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   &
697                  & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   &
698                  & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) )
699
700               ! negative part of the flux
701               zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   &
702                  & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   &
703                  & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) )
704
705               ! up & down beta terms
706               zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
707               zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
708               zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt
709            END DO
710         END DO
711      END DO
712      CALL lbc_lnk( zbetup, 'T', 1. )   ;   CALL lbc_lnk( zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign)
713
714      ! 3. monotonic flux in the i & j direction (paa & pbb)
715      ! ----------------------------------------
716      DO jk = 1, jpkm1
717         DO jj = 2, jpjm1
718            DO ji = fs_2, fs_jpim1   ! vector opt.
719               zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
720               zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
721               zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) )
722               paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
723
724               zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
725               zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
726               zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) )
727               pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
728
729      ! monotonic flux in the k direction, i.e. pcc
730      ! -------------------------------------------
731               za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
732               zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
733               zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
734               pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
735            END DO
736         END DO
737      END DO
738      CALL lbc_lnk( paa, 'U', -1. )   ;   CALL lbc_lnk( pbb, 'V', -1. )   ! lateral boundary condition (changed sign)
739      !
740      CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
741      !
742      IF( nn_timing == 1 )  CALL timing_stop('nonosc')
743      !
744   END SUBROUTINE nonosc
745
746
747   SUBROUTINE interp_4th_cpt( pt_in, pt_out )
748      !!----------------------------------------------------------------------
749      !!                  ***  ROUTINE interp_4th_cpt  ***
750      !!
751      !! **  Purpose :   Compute the interpolation of tracer at w-point
752      !!
753      !! **  Method  :   4th order compact interpolation
754      !!----------------------------------------------------------------------
755      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! now tracer fields
756      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! now tracer field interpolated at w-pts
757      !
758      INTEGER :: ji, jj, jk   ! dummy loop integers
759      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
760      !!----------------------------------------------------------------------
761     
762      DO jk = 3, jpkm1        !==  build the three diagonal matrix  ==!
763         DO jj = 1, jpj
764            DO ji = 1, jpi
765               zwd (ji,jj,jk) = 4._wp
766               zwi (ji,jj,jk) = 1._wp
767               zws (ji,jj,jk) = 1._wp
768               zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
769               !
770               IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom
771                  zwd (ji,jj,jk) = 1._wp
772                  zwi (ji,jj,jk) = 0._wp
773                  zws (ji,jj,jk) = 0._wp
774                  zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )   
775               ENDIF
776            END DO
777         END DO
778      END DO
779      !
780      jk=2                                            ! Switch to second order centered at top
781      DO jj=1,jpj
782         DO ji=1,jpi
783            zwd (ji,jj,jk) = 1._wp
784            zwi (ji,jj,jk) = 0._wp
785            zws (ji,jj,jk) = 0._wp
786            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
787         END DO
788      END DO   
789      !
790      !                       !==  tridiagonal solve  ==!
791      DO jj = 1, jpj                ! first recurrence
792         DO ji = 1, jpi
793            zwt(ji,jj,2) = zwd(ji,jj,2)
794         END DO
795      END DO
796      DO jk = 3, jpkm1
797         DO jj = 1, jpj
798            DO ji = 1, jpi
799               zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
800            END DO
801         END DO
802      END DO
803      !
804      DO jj = 1, jpj                ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
805         DO ji = 1, jpi
806            pt_out(ji,jj,2) = zwrm(ji,jj,2)
807         END DO
808      END DO
809      DO jk = 3, jpkm1
810         DO jj = 1, jpj
811            DO ji = 1, jpi
812               pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)             
813            END DO
814         END DO
815      END DO
816
817      DO jj = 1, jpj                ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
818         DO ji = 1, jpi
819            pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
820         END DO
821      END DO
822      DO jk = jpk-2, 2, -1
823         DO jj = 1, jpj
824            DO ji = 1, jpi
825               pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
826            END DO
827         END DO
828      END DO
829      !   
830   END SUBROUTINE interp_4th_cpt
831   
832   !!======================================================================
833END MODULE traadv_fct
Note: See TracBrowser for help on using the repository browser.