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_qck.F90 in NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_qck.F90 @ 13516

Last change on this file since 13516 was 13516, checked in by hadcv, 4 years ago

Tiling for tra_adv
  • Property svn:keywords set to Id
File size: 25.2 KB
Line 
1MODULE traadv_qck
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_qck  ***
4   !! Ocean tracers:  horizontal & vertical advective trend
5   !!==============================================================================
6   !! History :  3.0  !  2008-07  (G. Reffray)  Original code
7   !!            3.3  !  2010-05  (C.Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   tra_adv_qck    : update the tracer trend with the horizontal advection
12   !!                    trends using a 3rd order finite difference scheme
13   !!   tra_adv_qck_i  : apply QUICK scheme in i-direction
14   !!   tra_adv_qck_j  : apply QUICK scheme in j-direction
15   !!   tra_adv_cen2_k : 2nd centered scheme for the vertical advection
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and active tracers
18   USE dom_oce         ! ocean space and time domain
19   ! TEMP: This change not necessary after trd_tra is tiled
20   USE domain, ONLY : dom_tile
21   USE trc_oce         ! share passive tracers/Ocean variables
22   USE trd_oce         ! trends: ocean variables
23   USE trdtra          ! trends manager: tracers
24   USE diaptr          ! poleward transport diagnostics
25   USE iom
26   !
27   USE in_out_manager  ! I/O manager
28   USE lib_mpp         ! distribued memory computing
29   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
30   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   tra_adv_qck   ! routine called by step.F90
36
37   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
38
39   LOGICAL  ::   l_trd   ! flag to compute trends
40   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
41
42
43   !! * Substitutions
44#  include "do_loop_substitute.h90"
45#  include "domzgr_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
48   !! $Id$
49   !! Software governed by the CeCILL license (see ./LICENSE)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs )
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_adv_qck  ***
56      !!
57      !! ** Purpose :   Compute the now trend due to the advection of tracers
58      !!      and add it to the general trend of passive tracer equations.
59      !!
60      !! ** Method :   The advection is evaluated by a third order scheme
61      !!             For a positive velocity u :              u(i)>0
62      !!                                          |--FU--|--FC--|--FD--|------|
63      !!                                             i-1    i      i+1   i+2
64      !!
65      !!             For a negative velocity u :              u(i)<0
66      !!                                          |------|--FD--|--FC--|--FU--|
67      !!                                             i-1    i      i+1   i+2
68      !!             where  FU is the second upwind point
69      !!                    FD is the first douwning point
70      !!                    FC is the central point (or the first upwind point)
71      !!
72      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
73      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
74      !!
75      !!         dt = 2*rdtra and the scalar values are tb and sb
76      !!
77      !!       On the vertical, the simple centered scheme used pt(:,:,:,:,Kmm)
78      !!
79      !!               The fluxes are bounded by the ULTIMATE limiter to
80      !!             guarantee the monotonicity of the solution and to
81      !!            prevent the appearance of spurious numerical oscillations
82      !!
83      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends
84      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
85      !!             - poleward advective heat and salt transport (ln_diaptr=T)
86      !!
87      !! ** Reference : Leonard (1979, 1991)
88      !!----------------------------------------------------------------------
89      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index
90      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
91      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index
92      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
93      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers
94      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step
95      ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
96      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components
97      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation
98      !!----------------------------------------------------------------------
99      !
100      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
101         IF( kt == kit000 )  THEN
102            IF(lwp) WRITE(numout,*)
103            IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
104            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
105            IF(lwp) WRITE(numout,*)
106         ENDIF
107         !
108         l_trd = .FALSE.
109         l_ptr = .FALSE.
110         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
111         IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
112      ENDIF
113      !
114      !
115      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
116      CALL tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs ) 
117      CALL tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs ) 
118
119      !        ! vertical fluxes are computed with the 2nd order centered scheme
120      CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
121      !
122   END SUBROUTINE tra_adv_qck
123
124
125   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pU, Kbb, Kmm, pt, kjpt, Krhs )
126      !!----------------------------------------------------------------------
127      !!
128      !!----------------------------------------------------------------------
129      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
130      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
131      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
132      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
133      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
134      ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
135      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components
136      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
137      !!
138      ! TEMP: This change not necessary after trd_tra is tiled
139      INTEGER  ::   itile
140      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
141      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
142      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk) ::   zwx, zfu, zfc, zfd
143      ! TEMP: This change not necessary after trd_tra is tiled
144      REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   ztrdx
145      !----------------------------------------------------------------------
146      ! TEMP: This change not necessary after trd_tra is tiled
147      itile = ntile
148      !
149      ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
150      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
151         IF( kt == nit000 .AND. l_trd ) THEN
152            ALLOCATE( ztrdx(jpi,jpj,jpk) )
153         ENDIF
154      ENDIF
155      !                                                          ! ===========
156      DO jn = 1, kjpt                                            ! tracer loop
157         !                                                       ! ===========
158         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
159         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
160         !
161!!gm why not using a SHIFT instruction...
162         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
163            zfc(ji,jj,jk) = pt(ji-1,jj,jk,jn,Kbb)        ! Upstream   in the x-direction for the tracer
164            zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer
165         END_3D
166         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
167         
168         !
169         ! Horizontal advective fluxes
170         ! ---------------------------
171         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
172            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
173            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
174         END_3D
175         !
176         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
177            zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
178            zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm)
179            zwx(ji,jj,jk)  = ABS( pU(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
180            zfc(ji,jj,jk)  = zdir * pt(ji  ,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji+1,jj,jk,jn,Kbb)  ! FC in the x-direction for T
181            zfd(ji,jj,jk)  = zdir * pt(ji+1,jj,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji  ,jj,jk,jn,Kbb)  ! FD in the x-direction for T
182         END_3D
183         !--- Lateral boundary conditions
184         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp )
185
186         !--- QUICKEST scheme
187         CALL quickest( zfu, zfd, zfc, zwx )
188         !
189         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
190         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
191            zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
192         END_3D
193         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )      ! Lateral boundary conditions
194
195         !
196         ! Tracer flux on the x-direction
197         DO jk = 1, jpkm1 
198            !
199            DO_2D( 0, 0, 0, 0 )
200               zdir = 0.5 + SIGN( 0.5_wp, pU(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
201               !--- If the second ustream point is a land point
202               !--- the flux is computed by the 1st order UPWIND scheme
203               zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
204               zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
205               zwx(ji,jj,jk) = zwx(ji,jj,jk) * pU(ji,jj,jk)
206            END_2D
207         END DO
208         !
209         CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions
210         !
211         ! Computation of the trend
212         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
213            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
214            ! horizontal advective trends
215            ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
216            !--- add it to the general tracer trends
217            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
218         END_3D
219         !                                 ! trend diagnostics
220         ! TEMP: These changes not necessary after trd_tra is tiled
221         IF( l_trd )  THEN
222            DO_3D( 1, 0, 1, 0, 1, jpk )
223               ztrdx(ji,jj,jk) = zwx(ji,jj,jk)
224            END_3D
225
226            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
227               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
228
229               ! TODO: TO BE TILED- trd_tra
230               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) )
231
232               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = itile )     ! Revert to tile domain
233            ENDIF
234         ENDIF
235         !
236      END DO
237      !
238   END SUBROUTINE tra_adv_qck_i
239
240
241   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pV, Kbb, Kmm, pt, kjpt, Krhs )
242      !!----------------------------------------------------------------------
243      !!
244      !!----------------------------------------------------------------------
245      INTEGER                                  , INTENT(in   ) ::   kt         ! ocean time-step index
246      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
247      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
248      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers
249      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step
250      ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
251      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components
252      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
253      !!
254      ! TEMP: This change not necessary after trd_tra is tiled
255      INTEGER  ::   itile
256      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
257      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
258      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
259      ! TEMP: This change not necessary after trd_tra is tiled
260      REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   ztrdy
261      !----------------------------------------------------------------------
262      ! TEMP: This change not necessary after trd_tra is tiled
263      itile = ntile
264      !
265      ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
266      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
267         IF( kt == nit000 .AND. l_trd ) THEN
268            ALLOCATE( ztrdy(jpi,jpj,jpk) )
269         ENDIF
270      ENDIF
271      !                                                          ! ===========
272      DO jn = 1, kjpt                                            ! tracer loop
273         !                                                       ! ===========
274         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
275         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
276         !                                                 
277         DO jk = 1, jpkm1                               
278            !                                             
279            !--- Computation of the ustream and downstream value of the tracer and the mask
280            DO_2D( 0, 0, 0, 0 )
281               ! Upstream in the x-direction for the tracer
282               zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb)
283               ! Downstream in the x-direction for the tracer
284               zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb)
285            END_2D
286         END DO
287         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions
288
289         
290         !
291         ! Horizontal advective fluxes
292         ! ---------------------------
293         !
294         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
295            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
296            zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
297         END_3D
298         !
299         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
300            zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
301            zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
302            zwy(ji,jj,jk)  = ABS( pV(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
303            zfc(ji,jj,jk)  = zdir * pt(ji,jj  ,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj+1,jk,jn,Kbb)  ! FC in the x-direction for T
304            zfd(ji,jj,jk)  = zdir * pt(ji,jj+1,jk,jn,Kbb) + ( 1. - zdir ) * pt(ji,jj  ,jk,jn,Kbb)  ! FD in the x-direction for T
305         END_3D
306
307         !--- Lateral boundary conditions
308         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )
309
310         !--- QUICKEST scheme
311         CALL quickest( zfu, zfd, zfc, zwy )
312         !
313         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
314         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
315            zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
316         END_3D
317         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions
318         !
319         ! Tracer flux on the x-direction
320         DO jk = 1, jpkm1 
321            !
322            DO_2D( 0, 0, 0, 0 )
323               zdir = 0.5 + SIGN( 0.5_wp, pV(ji,jj,jk) )   ! if pU > 0 : zdir = 1 otherwise zdir = 0
324               !--- If the second ustream point is a land point
325               !--- the flux is computed by the 1st order UPWIND scheme
326               zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
327               zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
328               zwy(ji,jj,jk) = zwy(ji,jj,jk) * pV(ji,jj,jk)
329            END_2D
330         END DO
331         !
332         CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions
333         !
334         ! Computation of the trend
335         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
336            zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
337            ! horizontal advective trends
338            ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
339            !--- add it to the general tracer trends
340            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra
341         END_3D
342         !                                 ! trend diagnostics
343         ! TEMP: These changes not necessary after trd_tra is tiled
344         IF( l_trd )  THEN
345            DO_3D( 1, 0, 1, 0, 1, jpk )
346               ztrdy(ji,jj,jk) = zwy(ji,jj,jk)
347            END_3D
348
349            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
350               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
351
352               ! TODO: TO BE TILED- trd_tra
353               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) )
354
355               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = itile )     ! Revert to tile domain
356            ENDIF
357         ENDIF
358         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
359         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
360         !
361      END DO
362      !
363   END SUBROUTINE tra_adv_qck_j
364
365
366   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs )
367      !!----------------------------------------------------------------------
368      !!
369      !!----------------------------------------------------------------------
370      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
371      INTEGER                                  , INTENT(in   ) ::   Kmm, Krhs  ! ocean time level indices
372      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
373      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
374      ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
375      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity
376      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation
377      !
378      ! TEMP: This change not necessary after trd_tra is tiled
379      INTEGER  ::   itile
380      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
381      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk) ::   zwz   ! 3D workspace
382      ! TEMP: This change not necessary after trd_tra is tiled
383      REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   ztrdz
384      !!----------------------------------------------------------------------
385      ! TEMP: This change not necessary after trd_tra is tiled
386      itile = ntile
387      !
388      ! TEMP: This can be ST_2D(nn_hls) after trd_tra is tiled
389      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
390         IF( kt == nit000 .AND. l_trd ) THEN
391            ALLOCATE( ztrdz(jpi,jpj,jpk) )
392         ENDIF
393      ENDIF
394
395      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
396      zwz(:,:,jpk) = 0._wp
397      !
398      !                                                          ! ===========
399      DO jn = 1, kjpt                                            ! tracer loop
400         !                                                       ! ===========
401         !
402         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
403            zwz(ji,jj,jk) = 0.5 * pW(ji,jj,jk) * ( pt(ji,jj,jk-1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk)
404         END_3D
405         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
406            ! TODO: NOT TESTED- requires isf
407            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
408               DO_2D( 1, 1, 1, 1 )
409                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kmm)   ! linear free surface
410               END_2D
411            ELSE                                   ! no ocean cavities (only ocean surface)
412               DO_2D( 1, 1, 1, 1 )
413                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
414               END_2D
415            ENDIF
416         ENDIF
417         !
418         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
419            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
420               &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
421         END_3D
422         !                                 ! Send trends for diagnostic
423         ! TEMP: These changes not necessary after trd_tra is tiled
424         IF( l_trd )  THEN
425            DO_3D( 0, 0, 0, 0, 1, jpk )
426               ztrdz(ji,jj,jk) = zwz(ji,jj,jk)
427            END_3D
428
429            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
430               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )         ! Use full domain
431
432               ! TODO: TO BE TILED- trd_tra
433               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) )
434
435               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = itile )     ! Revert to tile domain
436            ENDIF
437         ENDIF
438         !
439      END DO
440      !
441   END SUBROUTINE tra_adv_cen2_k
442
443
444   SUBROUTINE quickest( pfu, pfd, pfc, puc )
445      !!----------------------------------------------------------------------
446      !!
447      !! ** Purpose :  Computation of advective flux with Quickest scheme
448      !!
449      !! ** Method :   
450      !!----------------------------------------------------------------------
451      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk), INTENT(in   ) ::   pfu   ! second upwind point
452      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk), INTENT(in   ) ::   pfd   ! first douwning point
453      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
454      REAL(wp), DIMENSION(ST_2D(nn_hls),jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
455      !!
456      INTEGER  ::  ji, jj, jk               ! dummy loop indices
457      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
458      REAL(wp) ::  zc, zcurv, zfho          !   -      -
459      !----------------------------------------------------------------------
460      !
461      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
462         zc     = puc(ji,jj,jk)                         ! Courant number
463         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
464         zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
465         zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
466         zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
467         zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
468         !
469         zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
470         zcoef2 = ABS( zcoef1 )
471         zcoef3 = ABS( zcurv )
472         IF( zcoef3 >= zcoef2 ) THEN
473            zfho = pfc(ji,jj,jk) 
474         ELSE
475            zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
476            IF( zcoef1 >= 0. ) THEN
477               zfho = MAX( pfc(ji,jj,jk), zfho ) 
478               zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
479            ELSE
480               zfho = MIN( pfc(ji,jj,jk), zfho ) 
481               zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
482            ENDIF
483         ENDIF
484         puc(ji,jj,jk) = zfho
485      END_3D
486      !
487   END SUBROUTINE quickest
488
489   !!======================================================================
490END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.