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 branches/dev_001_GM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 786

Last change on this file since 786 was 786, checked in by gm, 16 years ago

dev_001_GM - merge TRC-TRA on OPA only, trabbl & zpshde not done and trdmld not OK - compilation OK
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.2 KB
Line 
1MODULE traadv_qck
2   !!======================================================================
3   !!                       ***  MODULE  traadv_qck  ***
4   !! Ocean tracers:  horizontal & vertical advective trend using QUICKEST scheme
5   !!======================================================================
6   !! History :  1.0  !  06-09  (G. Reffray)  Original code
7   !!            2.4  !  08-01  (G. Madec)  Merge TRA-TRC
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   !!                 The vertical advection scheme is the 2nd centered scheme ?????
14   !!----------------------------------------------------------------------
15   USE dom_oce         ! ocean space and time domain
16   USE trdmod          ! ocean active tracers trends
17   USE trdmod_oce      ! ocean variables trends
18   USE flxrnf          !
19   USE ocfzpt          !
20   USE lib_mpp
21   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
22   USE in_out_manager  ! I/O manager
23   USE diaptr          ! poleward transport diagnostics
24   USE dynspg_oce      !
25   USE prtctl          ! Print control
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC tra_adv_qck    ! routine called by traadv.F90
31
32   REAL(wp), DIMENSION(jpi,jpj)     ::   zbtr2
33   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   sl
34   REAL(wp) ::   cst1, cst2, dt, coef1                  ! temporary scalars
35   INTEGER  ::   ji, jj, jk                             ! dummy loop indices
36
37   !! * Substitutions
38#  include "domzgr_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)
42   !! $Id:$
43   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
44   !!----------------------------------------------------------------------
45
46CONTAINS
47
48   SUBROUTINE tra_adv_qck( kt, cdtype, ktra, pun, pvn, pwn,   &
49      &                                      ptb, ptn, pta )
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE tra_adv_qck  ***
52      !!
53      !! ** Purpose :   Compute the now trend due to the advection of tracers
54      !!      and add it to the general trend of passive tracer equations.
55      !!
56      !! ** Method :   The advection is evaluated by a third order scheme
57      !!               For a positive velocity u :
58      !!
59      !!                  i-1    i      i+1   i+2
60      !!
61      !!               |--FU--|--FC--|--FD--|------|
62      !!                           u(i)>0
63      !!
64      !!               For a negative velocity u :
65      !!
66      !!               |------|--FD--|--FC--|--FU--|
67      !!                           u(i)<0
68      !!
69      !!                FU is the second upwind point
70      !!                FD is the first douwning point
71      !!                FC is the central point (or the first upwind point)
72      !!
73      !!      Flux(i) = u(i) * {0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i)Â?)/6)(FU+FD-2FC)}
74      !!                with C(i)=|u(i)|dx(i)/dt (Courant number)
75      !!
76      !!         dt = 2*rdtra and the scalar values are tb and sb
77      !!
78      !!       On the vertical, the simple centered scheme used tn and sn
79      !!
80      !!               The fluxes are bounded by the ULTIMATE limiter
81      !!               to guarantee the monotonicity of the solution and to
82      !!            prevent the appearance of spurious numerical oscillations
83      !!
84      !! ** Action : - update (ta,sa) with the now advective tracer trends
85      !!             - save the trends in (ttrdh,strdh) ('key_trdtra')
86      !!
87      !! References : Leonard (1979, 1991)
88      !!----------------------------------------------------------------------
89      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index
90      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator)
91      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index
92      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! 3 ocean velocity components
93      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn        ! before and now tracer fields
94      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend
95      !!
96      REAL(wp) :: z2                                   ! temporary scalar
97      !!----------------------------------------------------------------------
98
99      IF( kt == nit000 ) THEN
100         IF(lwp) WRITE(numout,*)
101         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3st order quickest advection scheme'
102         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
103
104         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
105         cst1 = 1./12.
106         cst2 = 2./3.
107      ENDIF
108
109      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.
110      ELSE                                        ;    z2 = 2.
111      ENDIF
112
113
114      ! I. Slope estimation at the T-point for the limiter ULTIMATE
115      ! SL = Sum(1/C_out) with C_out : Courant number for the outflows
116      !---------------------------------------------------------------
117
118      sl(:,:,:) = 100.
119
120                                                      ! ===============
121      DO jk = 1, jpkm1                                ! Horizontal slab
122        !                                             ! ===============
123        dt  = z2 * rdttra(jk)
124        DO jj = 2, jpjm1
125          DO ji = 2, fs_jpim1   ! vector opt.
126            coef1 = 1.e-12
127            IF (pun(ji-1,jj  ,jk  ).LT.0.) coef1 = coef1 + ABS(pun(ji-1,jj  ,jk  ))*dt/e1t(ji,jj)
128            IF (pun(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pun(ji  ,jj  ,jk  ))*dt/e1t(ji,jj)
129            IF (pvn(ji  ,jj-1,jk  ).LT.0.) coef1 = coef1 + ABS(pvn(ji  ,jj-1,jk  ))*dt/e2t(ji,jj)
130            IF (pvn(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pvn(ji  ,jj  ,jk  ))*dt/e2t(ji,jj)
131            IF (pwn(ji  ,jj  ,jk+1).LT.0.) coef1 = coef1 + ABS(pwn(ji  ,jj  ,jk+1))*dt/(fse3t(ji,jj,jk))
132            IF (pwn(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pwn(ji  ,jj  ,jk  ))*dt/(fse3t(ji,jj,jk))
133            sl(ji,jj,jk) = 1./coef1
134            sl(ji,jj,jk) = MIN(100.,sl(ji,jj,jk))
135            sl(ji,jj,jk) = MAX(1.  ,sl(ji,jj,jk))
136          ENDDO
137        ENDDO
138      ENDDO
139      !--- Lateral boundary conditions on the limiter slope
140      CALL lbc_lnk(   sl(:,:,:), 'T', 1. )
141
142
143      ! II. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
144      !---------------------------------------------------------------------------
145
146      CALL tra_adv_qck_hor( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 )
147
148      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - had: ', mask1=tmask, clinfo3=cdtype )
149
150
151      ! III. The vertical fluxes are computed with the 2nd order centered scheme
152      !-------------------------------------------------------------------------
153
154      CALL tra_adv_qck_ver( pwn, ptn , pta, z2 )
155
156      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - zad: ', mask1=tmask, clinfo3=cdtype )
157      !
158   END SUBROUTINE tra_adv_qck
159
160
161   SUBROUTINE tra_adv_qck_hor ( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 )
162      !!----------------------------------------------------------------------
163      !!
164      !!----------------------------------------------------------------------
165      INTEGER         , INTENT(in   )                         ::   kt          ! ocean time-step index
166      CHARACTER(len=3), INTENT(in   )                         ::   cdtype      ! =TRA or TRC (tracer indicator)
167      INTEGER         , INTENT(in   )                         ::   ktra        ! tracer index
168      REAL(wp)        , INTENT(in   )                         ::   z2          ! ???
169      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! horizontal effective velocity
170      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb         ! before tracer field
171      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta         ! tracer trend
172
173      REAL(wp) ::   za, zbtr, e1, e2, c, dir, fu, fc, fd   ! temporary scalars
174      REAL(wp) ::   coef2, coef3, fho, mask, dx
175      REAL(wp), DIMENSION(jpi,jpj) ::  zee
176      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmask, zlap, dwst, lim 
177      !!----------------------------------------------------------------------
178
179      !----------------------------------------------------------------------
180      ! 0. Initialization (should ot be needed on the whole array ???)
181      !----------------------------------------------------------------------
182      zmask(:,:,:)= 0.e0                                               
183      zlap (:,:,:)= 0.e0                                                 
184      dwst (:,:,:)= 0.e0                                                 
185      lim  (:,:,:)= 0.e0     
186                                           
187      !----------------------------------------------------------------------
188      ! I. Part 1 : x-direction
189      !----------------------------------------------------------------------
190
191                                                       ! ===============
192      DO jk = 1, jpkm1                                 ! Horizontal slab
193         !                                             ! ===============
194         ! Initialization of metric arrays (for z- or s-coordinates)
195         ! ---------------------------------------------------------
196         DO jj = 1, jpjm1
197            DO ji = 1, fs_jpim1   ! vector opt.
198#if defined key_zco
199               ! z-coordinates, no vertical scale factors
200               zee(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk)
201#else
202               ! vertical scale factor are used
203               zee(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
204#endif
205            END DO
206         END DO
207
208         ! Laplacian of tracers (at before time step)
209         ! ------------------------------------------
210         !--- First derivative (gradient)
211         DO jj = 1, jpjm1
212            DO ji = 1, fs_jpim1   ! vector opt.
213               zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) )
214            END DO
215         END DO
216         DO jj = 2, jpjm1
217            DO ji = fs_2, fs_jpim1   ! vector opt.
218#if defined key_zco
219               zee(ji,jj) = e1t(ji,jj) /  e2t(ji,jj)
220#else
221               zee(ji,jj) = e1t(ji,jj) / (e2t(ji,jj) * fse3t(ji,jj,jk))
222#endif
223               zlap(ji,jj,jk) = zee(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji-1,jj,jk)  )
224            END DO
225         END DO
226         !--- Function lim=FU+SL*(FC-FU) used by the limiter
227         !--- Computation of the ustream and downstream lim at the T-points
228!!gm bug : fs_2  instead of 2 ...
229!!gm  a lot of optimisation to be done in this routine....
230         DO jj = 2, jpjm1
231            DO ji = 2, fs_jpim1   ! vector opt.
232               ! Upstream in the x-direction for the tracer
233               zmask(ji,jj,jk) =ptb(ji-1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji-1,jj,jk))
234               ! Downstream in the x-direction for the tracer
235               dwst (ji,jj,jk) =ptb(ji+1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji+1,jj,jk))
236            ENDDO
237         ENDDO
238      END DO
239      !--- Lateral boundary conditions on the laplacian (unchanged sgn)
240      CALL lbc_lnk(   zlap(:,:,:), 'T', 1. ) 
241      !--- Lateral boundary conditions for the lim function
242      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  dwst(:,:,:), 'T', 1. )
243                                                       ! ===============
244      DO jk = 1, jpkm1                                 ! Horizontal slab
245         !                                             ! ===============
246         ! --- lim at the U-points in function of the direction of the flow
247         ! ----------------------------------------------------------------
248         DO jj = 1, jpjm1
249            DO ji = 1, fs_jpim1   ! vector opt.
250               dir = 0.5 + sign(0.5,pun(ji,jj,jk))     ! if pun>0 : diru = 1 otherwise diru = 0
251               lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji+1,jj,jk)
252               ! Mask at the T-points in the x-direction (mask=0 or mask=1)
253               zmask(ji,jj,jk)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2
254            END DO
255         END DO
256      END DO
257      !--- Lateral boundary conditions for the mask
258      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. ) 
259
260      ! Horizontal advective fluxes
261      ! ---------------------------
262                                                       ! ===============
263      DO jk = 1, jpkm1                                 ! Horizontal slab
264                                                       ! ===============
265         dt  = z2 * rdttra(jk)
266         !--- tracer flux at u and v-points
267         DO jj = 1, jpjm1
268            DO ji = 1, fs_jpim1   ! vector opt.
269#if defined key_zco
270               e2 = e2u(ji,jj)
271#else
272               e2 = e2u(ji,jj) * fse3u(ji,jj,jk)
273#endif
274               dir = 0.5 + sign(0.5,pun(ji,jj,jk))                ! if pun>0 : diru = 1 otherwise diru = 0
275
276               dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj)
277               c  = ABS(pun(ji,jj,jk))*dt/dx                      ! (0<cx<1 : Courant number on x-direction)
278
279               fu  = lim(ji,jj,jk)                                ! FU + sl(FC-FU) in the x-direction for T
280               fc  = dir*ptb(ji  ,jj,jk)+(1-dir)*ptb(ji+1,jj,jk)  ! FC in the x-direction for T
281               fd  = dir*ptb(ji+1,jj,jk)+(1-dir)*ptb(ji  ,jj,jk)  ! FD in the x-direction for T
282
283               !--- QUICKEST scheme
284               ! Temperature on the x-direction
285               coef1 = 0.5*(fc+fd)
286               coef2 = 0.5*c*(fd-fc)
287               coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji+1,jj,jk) )
288               fho = coef1-coef2-coef3
289               fho = bound(fu,fd,fc,fho)
290               !--- If the second ustream point is a land point
291               !--- the flux is computed by the 1st order UPWIND scheme
292               mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji+1,jj,jk)
293               fho = mask*fho + (1-mask)*fc
294               dwst(ji,jj,jk)=e2*pun(ji,jj,jk)*fho
295            END DO
296         END DO
297
298         !--- Tracer flux divergence at t-point added to the general trend
299         DO jj = 2, jpjm1
300            DO ji = fs_2, fs_jpim1   ! vector opt.
301               !--- horizontal advective trends
302#if defined key_zco
303               zbtr = zbtr2(ji,jj)
304#else
305               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)             
306#endif
307               za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji-1,jj  ,jk) )
308               !--- add it to the general tracer trends
309               pta(ji,jj,jk) = pta(ji,jj,jk) + za
310            END DO
311         END DO
312         !                                             ! ===============
313      END DO                                           !   End of slab
314      !                                                ! ===============
315
316      ! Save the horizontal advective trends for diagnostic
317      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, dwst, pun, ptb )
318
319
320      !----------------------------------------------------------------------
321      ! I. Part 2 : y-direction
322      !----------------------------------------------------------------------
323                                                       ! ==============
324      DO jk = 1, jpkm1                                 ! Horizontal slab
325         !                                             ! ===============
326         ! Initialization of metric arrays (for z- or s-coordinates)
327         ! ---------------------------------------------------------
328         DO jj = 1, jpjm1
329            DO ji = 1, fs_jpim1   ! vector opt.
330#if defined key_zco
331               ! z-coordinates, no vertical scale factors
332               zee(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk)
333#else
334               ! s-coordinates, vertical scale factor are used
335               zee(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
336#endif
337            END DO
338         END DO
339
340         ! Laplacian of tracers (at before time step)
341         ! ------------------------------------------
342         !--- First derivative (gradient)
343         DO jj = 1, jpjm1
344            DO ji = 1, fs_jpim1   ! vector opt.
345               zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) )
346            END DO
347         END DO
348         !--- Second derivative (divergence)
349         DO jj = 2, jpjm1
350            DO ji = fs_2, fs_jpim1   ! vector opt.
351#if defined key_zco
352               zee(ji,jj) = e2t(ji,jj) /  e1t(ji,jj)
353#else
354               zee(ji,jj) = e2t(ji,jj) / (e1t(ji,jj) * fse3t(ji,jj,jk))
355#endif
356               zlap(ji,jj,jk) = zee(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji,jj-1,jk)  )
357            END DO
358         END DO
359         !--- Function lim=FU+SL*(FC-FU) used by the limiter
360         !--- Computation of the ustream and downstream lim at the T-points
361         DO jj = 2, jpjm1
362            DO ji = 2, fs_jpim1   ! vector opt.
363               ! Upstream in the y-direction for the tracer
364               zmask(ji,jj,jk)=ptb(ji,jj-1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj-1,jk))
365               ! Downstream in the y-direction for the tracer
366               dwst (ji,jj,jk)=ptb(ji,jj+1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj+1,jk))
367            ENDDO
368         ENDDO
369      END DO
370      !--- Lateral boundary conditions on the laplacian (unchanged sgn)
371      CALL lbc_lnk(   zlap(:,:,:), 'T', 1. )
372      !--- Lateral boundary conditions for the lim function
373      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  dwst(:,:,:), 'T', 1. )
374
375      DO jk = 1, jpkm1                                 ! Horizontal slab
376         !                                             ! ===============
377         ! --- lim at the V-points in function of the direction of the flow
378         ! ----------------------------------------------------------------
379         DO jj = 1, jpjm1
380            DO ji = 1, fs_jpim1   ! vector opt.
381               dir = 0.5 + sign(0.5,pvn(ji,jj,jk))      ! if pvn>0 : dirv = 1 otherwise dirv = 0
382               lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji,jj+1,jk)
383               ! Mask at the T-points in the y-direction (mask=0 or mask=1)
384               zmask(ji,jj,jk)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2
385            END DO
386         END DO
387      END DO
388      !--- Lateral boundary conditions for the mask
389      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. ) 
390
391      ! Horizontal advective fluxes
392      ! -------------------------------
393                                                       ! ===============
394      DO jk = 1, jpkm1                                 ! Horizontal slab
395                                                       ! ===============
396         dt  = z2 * rdttra(jk)
397         !--- tracer flux at u and v-points
398         DO jj = 1, jpjm1
399            DO ji = 1, fs_jpim1   ! vector opt.
400#if defined key_zco
401               e1 = e1v(ji,jj)
402#else
403               e1 = e1v(ji,jj) * fse3v(ji,jj,jk)
404#endif
405               dir = 0.5 + sign(0.5,pvn(ji,jj,jk))                ! if pvn>0 : dirv = 1 otherwise dirv = 0
406
407               dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1)
408               c  = ABS(pvn(ji,jj,jk))*dt/dx                      ! (0<cy<1 : Courant number on y-direction)
409
410               fu  = lim(ji,jj,jk)                                ! FU + sl(FC-FU) in the y-direction for T
411               fc  = dir*ptb(ji,jj  ,jk)+(1-dir)*ptb(ji,jj+1,jk)  ! FC in the y-direction for T
412               fd  = dir*ptb(ji,jj+1,jk)+(1-dir)*ptb(ji,jj  ,jk)  ! FD in the y-direction for T
413
414               !--- QUICKEST scheme
415               ! Temperature on the y-direction
416               coef1 = 0.5*(fc+fd)
417               coef2 = 0.5*c*(fd-fc)
418               coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji,jj+1,jk) )
419               fho = coef1-coef2-coef3
420               fho = bound(fu,fd,fc,fho)
421               !--- If the second ustream point is a land point
422               !--- the flux is computed by the 1st order UPWIND scheme
423               mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji,jj+1,jk)
424               fho = mask*fho + (1-mask)*fc
425               dwst(ji,jj,jk)=e1*pvn(ji,jj,jk)*fho
426            END DO
427         END DO
428
429         !--- Tracer flux divergence at t-point added to the general trend
430         DO jj = 2, jpjm1
431            DO ji = fs_2, fs_jpim1   ! vector opt.
432               !--- horizontal advective trends
433#if defined key_zco
434               zbtr = zbtr2(ji,jj)
435#else
436               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
437#endif
438               za = - zbtr * (  dwst(ji,jj,jk) - dwst(ji  ,jj-1,jk) )
439               !--- add it to the general tracer trends
440               pta(ji,jj,jk) = pta(ji,jj,jk) + za
441            END DO
442         END DO
443         !                                             ! ===============
444      END DO                                           !   End of slab
445      !                                                ! ===============
446
447      ! Save the horizontal advective trends for diagnostic
448      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, dwst, pvn, ptb )
449
450      ! "zonal" mean advective heat and salt transport
451      IF( ln_diaptr .AND. cdtype == 'TRA' .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
452#if defined key_zco
453         DO jk = 1, jpkm1
454            DO jj = 2, jpjm1
455               DO ji = fs_2, fs_jpim1   ! vector opt.
456                  dwst(ji,jj,jk) = dwst(ji,jj,jk) * fse3v(ji,jj,jk)
457               END DO
458            END DO
459         END DO
460# endif
461         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( dwst(:,:,:) )
462         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( dwst(:,:,:) )
463      ENDIF
464      !
465   END SUBROUTINE tra_adv_qck_hor
466
467
468   SUBROUTINE tra_adv_qck_ver( pwn, ptn , pta, z2  )     
469      !!----------------------------------------------------------------------
470      !!
471      !!----------------------------------------------------------------------
472      REAL(wp), INTENT(in   )                         ::   z2
473      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn
474      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptn
475      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta
476      !!
477      REAL(wp) ::   za, ze3tr, dt, dir, fc, fd             ! temporary scalars
478
479      ! Vertical advection
480      ! ------------------
481
482      ! 1. Vertical advective fluxes
483      ! ----------------------------
484
485      sl(:,:,jpk) = 0.e0      !Bottom value : flux set to zero
486
487      ! Surface value
488      IF( lk_dynspg_rl .OR. lk_vvl ) THEN         ! rigid lid : flux set to zero
489         sl(:,:, 1 ) = 0.e0
490      ELSE                                        ! free surface-constant volume
491         sl(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)
492      ENDIF
493
494      ! Second order centered tracer flux at w-point
495      DO jk = 2, jpkm1
496         dt  = z2 *  rdttra(jk)
497         DO jj = 2, jpjm1
498            DO ji = fs_2, fs_jpim1   ! vector opt.
499               dir = 0.5 + sign(0.5,pwn(ji,jj,jk))                                         ! if pwn>0 : dirw = 1 otherwise dirw = 0
500               fc = dir*ptn(ji,jj,jk  )*fse3t(ji,jj,jk-1)+(1-dir)*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk  )   ! FC in the z-direction for T
501               fd = dir*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk  )+(1-dir)*ptn(ji,jj,jk  )*fse3t(ji,jj,jk-1)   ! FD in the z-direction for T
502               !--- Second order centered scheme
503               sl(ji,jj,jk)=pwn(ji,jj,jk)*(fc+fd)/(fse3t(ji,jj,jk-1)+fse3t(ji,jj,jk))
504            END DO
505         END DO
506      END DO
507
508      ! 2. Tracer flux divergence at t-point added to the general trend
509      ! ---------------------------------------------------------------
510      DO jk = 1, jpkm1
511         DO jj = 2, jpjm1
512            DO ji = fs_2, fs_jpim1   ! vector opt.
513               ze3tr = 1. / fse3t(ji,jj,jk)
514               ! vertical advective trends
515               za = - ze3tr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) )
516               ! add it to the general tracer trends
517               pta(ji,jj,jk) =  pta(ji,jj,jk) + za
518            END DO
519         END DO
520      END DO
521      !
522   END SUBROUTINE tra_adv_qck_ver
523
524
525   REAL FUNCTION bound(fu,fd,fc,fho)
526      REAL(wp) ::  fu, fd, fc, fho, fref1, fref2
527      fref1 = fu
528      fref2 = MAX(  MIN( fc , fd ), MIN( MAX( fc , fd ), fref1 )  )
529      bound = MAX(  MIN( fho, fc ), MIN( MAX( fho, fc ), fref2 )  )
530   END FUNCTION
531
532   !!======================================================================
533END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.