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/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/OCE/TRA – NEMO

source: NEMO/branches/UKMO/dev_r12866_HPC-02_Daley_Tiling_trial_extra_halo/src/OCE/TRA/traadv_qck.F90 @ 13409

Last change on this file since 13409 was 13409, checked in by hadcv, 3 years ago

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