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

source: branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 @ 6995

Last change on this file since 6995 was 6995, checked in by acc, 7 years ago

Branch dev_r6393_NOC_WAD. Corrected reproducibility error in wet_dry.F90

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