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

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 9 years ago

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