source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 9 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

  • Property svn:keywords set to Id
File size: 29.5 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   USE trdmod_oce      ! ocean space and time domain
20   USE trdtra          ! ocean tracers trends
21   USE trabbl          ! advective term in the BBL
22   USE lib_mpp         ! distribued memory computing
23   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
24   USE dynspg_oce      ! surface pressure gradient variables
25   USE in_out_manager  ! I/O manager
26   USE diaptr          ! poleward transport diagnostics
27   USE trc_oce         ! share passive tracers/Ocean variables
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC   tra_adv_qck   ! routine called by step.F90
33
34   LOGICAL  :: l_trd           ! flag to compute trends
35   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
36
37   !! * Control permutation of array indices
38#  include "oce_ftrans.h90"
39#  include "dom_oce_ftrans.h90"
40#  include "trc_oce_ftrans.h90"
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE tra_adv_qck ( kt, cdtype, p2dt, pun, pvn, pwn,      &
53      &                                       ptb, ptn, pta, kjpt )
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 ptn
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 (pta) with the now advective tracer trends
84      !!             - save the trends
85      !!
86      !! ** Reference : Leonard (1979, 1991)
87      !!----------------------------------------------------------------------
88      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
89      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
90      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
91      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
92
93      !! DCSE_NEMO: This style defeats ftrans
94!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
95!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
96!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
97
98!FTRANS pun pvn pwn :I :I :z
99!FTRANS ptb ptn :I :I :z :
100!FTRANS pta :I :I :z :
101      REAL(wp), INTENT(in   ) ::   pun(jpi,jpj,jpk)        ! ocean velocity component (u)
102      REAL(wp), INTENT(in   ) ::   pvn(jpi,jpj,jpk)        ! ocean velocity component (v)
103      REAL(wp), INTENT(in   ) ::   pwn(jpi,jpj,jpk)        ! ocean velocity component (w)
104      REAL(wp), INTENT(in   ) ::   ptb(jpi,jpj,jpk,kjpt)   ! tracer fields (before)
105      REAL(wp), INTENT(in   ) ::   ptn(jpi,jpj,jpk,kjpt)   ! tracer fields (now)
106      REAL(wp), INTENT(inout) ::   pta(jpi,jpj,jpk,kjpt)   ! tracer trend
107
108      !!----------------------------------------------------------------------
109
110      IF( kt == nit000 )  THEN
111         IF(lwp) WRITE(numout,*)
112         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
113         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
114         IF(lwp) WRITE(numout,*)
115         !
116         l_trd = .FALSE.
117         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
118      ENDIF
119
120      ! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
121      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
122      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
123
124      ! II. The vertical fluxes are computed with the 2nd order centered scheme
125      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
126      !
127
128      !! * Reset control of array index permutation
129!FTRANS CLEAR
130#  include "oce_ftrans.h90"
131#  include "dom_oce_ftrans.h90"
132#  include "trc_oce_ftrans.h90"
133
134   END SUBROUTINE tra_adv_qck
135
136
137   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
138      &                                        ptb, ptn, pta, kjpt   )
139      !!----------------------------------------------------------------------
140      !!
141      !!----------------------------------------------------------------------
142      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
143      USE oce     , ONLY:   zwx => ua       ! ua used as workspace
144      USE wrk_nemo, ONLY:   zfu => wrk_3d_1 , zfc => wrk_3d_2, zfd => wrk_3d_3   ! 3D workspace
145
146      !! DCSE_NEMO: need additional directives for renamed module variables
147!FTRANS zwx zfu zfc zfd :I :I :z
148
149      !
150      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
151      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
152      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
153      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
154
155      !! DCSE_NEMO: This style defeats ftrans
156!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-velocity components
157!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
158!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
159
160!FTRANS pun :I :I :z
161!FTRANS ptb ptn pta :I :I :z :
162      REAL(wp), INTENT(in   ) ::   pun(jpi,jpj,jpk)         ! i-velocity component
163      REAL(wp), INTENT(in   ) ::   ptb(jpi,jpj,jpk,kjpt)    ! tracer field (before)
164      REAL(wp), INTENT(in   ) ::   ptn(jpi,jpj,jpk,kjpt)    ! tracer field (now)
165      REAL(wp), INTENT(inout) ::   pta(jpi,jpj,jpk,kjpt)    ! tracer trend
166
167      !!
168      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
169      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
170      !----------------------------------------------------------------------
171      !
172      IF( wrk_in_use(3, 1,2,3) ) THEN
173         CALL ctl_stop('tra_adv_qck_i: requested workspace arrays unavailable')   ;   RETURN
174      ENDIF
175      !                                                          ! ===========
176      DO jn = 1, kjpt                                            ! tracer loop
177         !                                                       ! ===========
178         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
179         zfd(:,:,:) = 0.0     ;   zwx(:,:,:) = 0.0     
180         !                                                 
181#if defined key_z_first
182         !--- Computation of the upstream and downstream value of the tracer and the mask
183         DO jj = 2, jpjm1
184            DO ji = 2, jpim1
185               DO jk = 1, jpkm1                               
186#else
187         DO jk = 1, jpkm1                               
188            !                                             
189            !--- Computation of the upstream and downstream value of the tracer and the mask
190            DO jj = 2, jpjm1
191               DO ji = fs_2, fs_jpim1   ! vector opt.
192#endif
193                  ! Upstream in the x-direction for the tracer
194                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
195                  ! Downstream in the x-direction for the tracer
196                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
197               END DO
198            END DO
199         END DO
200         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
201         
202         !
203         ! Horizontal advective fluxes
204         ! ---------------------------
205         !
206#if defined key_z_first
207         DO jj = 2, jpjm1
208            DO ji = 2, jpim1
209               DO jk = 1, jpkm1                             
210#else
211         DO jk = 1, jpkm1                             
212            DO jj = 2, jpjm1
213               DO ji = fs_2, fs_jpim1   ! vector opt.         
214#endif
215                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
216                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
217               END DO
218            END DO
219         END DO
220         !
221#if defined key_z_first
222         DO jj = 2, jpjm1
223            DO ji = 2, jpim1
224               DO jk = 1, jpkm1 
225                  zdt =  p2dt(jk)
226#else
227         DO jk = 1, jpkm1 
228            zdt =  p2dt(jk)
229            DO jj = 2, jpjm1
230               DO ji = fs_2, fs_jpim1   ! vector opt.   
231#endif
232                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
233                  zdx  = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
234                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
235                  zfc(ji,jj,jk)  = zdir * ptb(ji  ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn)  ! FC in the x-direction for T
236                  zfd(ji,jj,jk)  = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji  ,jj,jk,jn)  ! FD in the x-direction for T
237               END DO
238            END DO
239         END DO 
240         !--- Lateral boundary conditions
241         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
242         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
243
244         !--- QUICKEST scheme
245         CALL quickest( zfu, zfd, zfc, zwx )
246         !
247         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
248#if defined key_z_first
249         DO jj = 2, jpjm1
250            DO ji = 2, jpim1
251               DO jk = 1, jpkm1 
252#else
253         DO jk = 1, jpkm1 
254            DO jj = 2, jpjm1
255               DO ji = fs_2, fs_jpim1   ! vector opt.               
256#endif
257                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
258               END DO
259            END DO
260         END DO
261         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
262
263         !
264         ! Tracer flux on the x-direction
265#if defined key_z_first
266         DO jj = 2, jpjm1
267            DO ji = 2, jpim1
268               DO jk = 1, jpkm1 
269#else
270         DO jk = 1, jpkm1 
271            DO jj = 2, jpjm1
272               DO ji = fs_2, fs_jpim1   ! vector opt.               
273#endif
274                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
275                  !--- If the second ustream point is a land point
276                  !--- the flux is computed by the 1st order UPWIND scheme
277                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
278                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
279                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
280               END DO
281            END DO
282#if defined key_z_first
283         END DO
284         ! Computation of the trend
285         DO jj = 2, jpjm1
286            DO ji = 2, jpim1
287               DO jk = 1, jpkm1
288#else
289            !
290            ! Computation of the trend
291            DO jj = 2, jpjm1
292               DO ji = fs_2, fs_jpim1   ! vector opt. 
293#endif
294                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
295                  ! horizontal advective trends
296                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
297                  !--- add it to the general tracer trends
298                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
299               END DO
300            END DO
301            !
302         END DO
303         !                                 ! trend diagnostics (contribution of upstream fluxes)
304         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
305         !
306      END DO
307      !
308      IF( wrk_not_released(3, 1,2,3) )   CALL ctl_stop('tra_adv_qck_i: failed to release workspace arrays')
309      !
310
311      !! * Reset control of array index permutation
312!FTRANS CLEAR
313#  include "oce_ftrans.h90"
314#  include "dom_oce_ftrans.h90"
315#  include "trc_oce_ftrans.h90"
316
317   END SUBROUTINE tra_adv_qck_i
318
319
320   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
321      &                                        ptb, ptn, pta, kjpt )
322      !!----------------------------------------------------------------------
323      !!
324      !!----------------------------------------------------------------------
325      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
326      USE oce     , ONLY:   zwy => ua       ! ua used as workspace
327      USE wrk_nemo, ONLY:   zfu => wrk_3d_1 , zfc => wrk_3d_2, zfd => wrk_3d_3   ! 3D workspace
328
329      !! DCSE_NEMO: need additional directives for renamed module variables
330!FTRANS zwy zfu zfc zfd :I :I :z
331
332      !
333      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
334      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
335      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
336      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
337
338      !! DCSE_NEMO: This style defeats ftrans
339!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
340!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
341!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
342
343!FTRANS pvn :I :I :z
344!FTRANS ptb ptn pta :I :I :z :
345      REAL(wp), INTENT(in   ) ::   pvn(jpi,jpj,jpk)         ! j-velocity component
346      REAL(wp), INTENT(in   ) ::   ptb(jpi,jpj,jpk,kjpt)    ! tracer field (before)
347      REAL(wp), INTENT(in   ) ::   ptn(jpi,jpj,jpk,kjpt)    ! tracer field (now)
348      REAL(wp), INTENT(inout) ::   pta(jpi,jpj,jpk,kjpt)    ! tracer trend
349
350      !!
351      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
352      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
353      !----------------------------------------------------------------------
354      !
355      IF(wrk_in_use(3, 1,2,3))THEN
356         CALL ctl_stop('tra_adv_qck_j: ERROR: requested workspace arrays unavailable')
357         RETURN
358      END IF
359      !                                                          ! ===========
360      DO jn = 1, kjpt                                            ! tracer loop
361         !                                                       ! ===========
362         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
363         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
364         !                                                 
365#if defined key_z_first
366         !--- Computation of the ustream and downstream value of the tracer and the mask
367         DO jj = 2, jpjm1
368            DO ji = 2, jpim1
369               DO jk = 1, jpkm1                               
370#else
371         DO jk = 1, jpkm1                               
372            !                                             
373            !--- Computation of the ustream and downstream value of the tracer and the mask
374            DO jj = 2, jpjm1
375               DO ji = fs_2, fs_jpim1   ! vector opt.
376#endif
377                  ! Upstream in the x-direction for the tracer
378                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
379                  ! Downstream in the x-direction for the tracer
380                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
381               END DO
382            END DO
383         END DO
384         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
385
386         
387         !
388         ! Horizontal advective fluxes
389         ! ---------------------------
390         !
391#if defined key_z_first
392         DO jj = 2, jpjm1
393            DO ji = 2, jpim1
394               DO jk = 1, jpkm1                             
395#else
396         DO jk = 1, jpkm1                             
397            DO jj = 2, jpjm1
398               DO ji = fs_2, fs_jpim1   ! vector opt.         
399#endif
400                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
401                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
402               END DO
403            END DO
404         END DO
405         !
406#if defined key_z_first
407         DO jj = 2, jpjm1
408            DO ji = 2, jpim1
409               DO jk = 1, jpkm1 
410                  zdt =  p2dt(jk)
411#else
412         DO jk = 1, jpkm1 
413            zdt =  p2dt(jk)
414            DO jj = 2, jpjm1
415               DO ji = fs_2, fs_jpim1   ! vector opt.   
416#endif
417                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
418                  zdx  = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
419                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
420                  zfc(ji,jj,jk)  = zdir * ptb(ji,jj  ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn)  ! FC in the x-direction for T
421                  zfd(ji,jj,jk)  = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj  ,jk,jn)  ! FD in the x-direction for T
422               END DO
423            END DO
424         END DO
425
426         !--- Lateral boundary conditions
427         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
428         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
429
430         !--- QUICKEST scheme
431         CALL quickest( zfu, zfd, zfc, zwy )
432         !
433         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
434#if defined key_z_first
435         DO jj = 2, jpjm1
436            DO ji = 2, jpim1
437               DO jk = 1, jpkm1 
438#else
439         DO jk = 1, jpkm1 
440            DO jj = 2, jpjm1
441               DO ji = fs_2, fs_jpim1   ! vector opt.               
442#endif
443                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
444               END DO
445            END DO
446         END DO
447         !--- Lateral boundary conditions
448         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
449         !
450         ! Tracer flux on the x-direction
451#if defined key_z_first
452         DO jj = 2, jpjm1
453            DO ji = 2, jpim1
454               DO jk = 1, jpkm1 
455#else
456         DO jk = 1, jpkm1 
457            !
458            DO jj = 2, jpjm1
459               DO ji = fs_2, fs_jpim1   ! vector opt.               
460#endif
461                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
462                  !--- If the second ustream point is a land point
463                  !--- the flux is computed by the 1st order UPWIND scheme
464                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
465                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
466                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
467               END DO
468            END DO
469#if defined key_z_first
470         END DO
471         ! Computation of the trend
472         DO jj = 2, jpjm1
473            DO ji = 2, jpim1
474               DO jk = 1, jpkm1
475#else
476            !
477            ! Computation of the trend
478            DO jj = 2, jpjm1
479               DO ji = fs_2, fs_jpim1   ! vector opt. 
480#endif
481                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
482                  ! horizontal advective trends
483                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
484                  !--- add it to the general tracer trends
485                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
486               END DO
487            END DO
488            !
489         END DO
490         !                                 ! trend diagnostics (contribution of upstream fluxes)
491         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
492         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
493         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
494           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
495           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
496         ENDIF
497         !
498      END DO
499      !
500      IF( wrk_not_released(3, 1,2,3) )   CALL ctl_stop('tra_adv_qck_j: failed to release workspace arrays')
501      !
502
503      !! * Reset control of array index permutation
504!FTRANS CLEAR
505#  include "oce_ftrans.h90"
506#  include "dom_oce_ftrans.h90"
507#  include "trc_oce_ftrans.h90"
508
509   END SUBROUTINE tra_adv_qck_j
510
511
512   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
513     &                                    ptn, pta, kjpt )
514      !!----------------------------------------------------------------------
515      !!
516      !!----------------------------------------------------------------------
517      USE oce, ONLY:   zwz => ua   ! ua used as workspace
518
519      !! DCSE_NEMO: need additional directives for renamed module variables
520!FTRANS zwz :I :I :z
521
522      !
523      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
524      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
525      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
526
527      !! DCSE_NEMO: This style defeats ftrans
528!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
529!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! tracer fields (now)
530!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
531
532!FTRANS pwn :I :I :z
533!FTRANS ptn pta :I :I :z :
534      REAL(wp), INTENT(in   ) ::   pwn(jpi,jpj,jpk)         ! vertical velocity
535      REAL(wp), INTENT(in   ) ::   ptn(jpi,jpj,jpk,kjpt)    ! tracer fields (now)
536      REAL(wp), INTENT(inout) ::   pta(jpi,jpj,jpk,kjpt)    ! tracer trend
537
538      !
539      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
540      REAL(wp) ::   zbtr , ztra      ! local scalars
541      !!----------------------------------------------------------------------
542
543      !                                                          ! ===========
544      DO jn = 1, kjpt                                            ! tracer loop
545         !                                                       ! ===========
546         ! 1. Bottom value : flux set to zero
547         zwz(:,:,jpk) = 0.e0             ! Bottom value : flux set to zero
548         !
549         !                                 ! Surface value
550         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! Variable volume : flux set to zero
551         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! Constant volume : advective flux through the surface
552         ENDIF
553         !
554#if defined key_z_first
555         DO jj = 2, jpjm1
556            DO ji = 2, jpim1
557               DO jk = 2, jpkm1            ! Interior point: second order centered tracer flux at w-point
558#else
559         DO jk = 2, jpkm1                  ! Interior point: second order centered tracer flux at w-point
560            DO jj = 2, jpjm1
561               DO ji = fs_2, fs_jpim1   ! vector opt.
562#endif
563                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
564               END DO
565            END DO
566         END DO
567         !
568#if defined key_z_first
569         DO jj = 2, jpjm1
570            DO ji = 2, jpim1
571               DO jk = 1, jpkm1    !==  Tracer flux divergence added to the general trend  ==!
572#else
573         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
574            DO jj = 2, jpjm1
575               DO ji = fs_2, fs_jpim1   ! vector opt.
576#endif
577                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
578                  ! k- vertical advective trends
579                  ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 
580                  ! added to the general tracer trends
581                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
582               END DO
583            END DO
584         END DO
585         !                                 ! Save the vertical advective trends for diagnostic
586         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
587         !
588      END DO
589      !
590
591      !! * Reset control of array index permutation
592!FTRANS CLEAR
593#  include "oce_ftrans.h90"
594#  include "dom_oce_ftrans.h90"
595#  include "trc_oce_ftrans.h90"
596
597   END SUBROUTINE tra_adv_cen2_k
598
599
600   SUBROUTINE quickest( pfu, pfd, pfc, puc )
601      !!----------------------------------------------------------------------
602      !!
603      !! ** Purpose :  Computation of advective flux with Quickest scheme
604      !!
605      !! ** Method :   
606      !!----------------------------------------------------------------------
607
608      !! DCSE_NEMO: This style defeats ftrans
609
610!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
611!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
612!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
613!     REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
614
615!FTRANS pfu pfd pfc puc :I :I :z
616      REAL(wp), INTENT(in   ) ::   pfu(jpi,jpj,jpk)   ! second upwind point
617      REAL(wp), INTENT(in   ) ::   pfd(jpi,jpj,jpk)   ! first douwning point
618      REAL(wp), INTENT(in   ) ::   pfc(jpi,jpj,jpk)   ! the central point (or the first upwind point)
619      REAL(wp), INTENT(inout) ::   puc(jpi,jpj,jpk)   ! input as Courant number ; output as flux
620
621      !!
622      INTEGER  ::  ji, jj, jk               ! dummy loop indices
623      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
624      REAL(wp) ::  zc, zcurv, zfho          !   -      -
625      !----------------------------------------------------------------------
626
627#if defined key_z_first
628      DO jj = 1, jpj
629         DO ji = 1, jpi
630            DO jk = 1, jpkm1
631#else
632      DO jk = 1, jpkm1
633         DO jj = 1, jpj
634            DO ji = 1, jpi
635#endif
636               zc     = puc(ji,jj,jk)                         ! Courant number
637               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
638               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
639               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
640               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
641               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
642               !
643               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
644               zcoef2 = ABS( zcoef1 )
645               zcoef3 = ABS( zcurv )
646               IF( zcoef3 >= zcoef2 ) THEN
647                  zfho = pfc(ji,jj,jk) 
648               ELSE
649                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
650                  IF( zcoef1 >= 0. ) THEN
651                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
652                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
653                  ELSE
654                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
655                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
656                  ENDIF
657               ENDIF
658               puc(ji,jj,jk) = zfho
659            END DO
660         END DO
661      END DO
662      !
663   END SUBROUTINE quickest
664
665   !!======================================================================
666END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.