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

source: branches/2015/dev_CMCC_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 @ 6051

Last change on this file since 6051 was 6051, checked in by lovato, 8 years ago

Merge branches/2015/dev_r5056_CMCC4_simplification (see ticket #1456)

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