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

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 5883

Last change on this file since 5883 was 5883, checked in by gm, 8 years ago

#1613: vvl by default: TRA/TRC remove optimization associated with linear free surface
  • Property svn:keywords set to Id
File size: 24.1 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 trc_oce         ! share passive tracers/Ocean variables
20   USE trd_oce         ! trends: ocean variables
21   USE trdtra          ! trends manager: tracers
22   USE dynspg_oce      ! surface pressure gradient variables
23   USE diaptr          ! poleward transport diagnostics
24   !
25   USE lib_mpp         ! distribued memory computing
26   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
27   USE in_out_manager  ! I/O manager
28   USE wrk_nemo        ! Memory Allocation
29   USE timing          ! Timing
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   LOGICAL  :: l_trd           ! flag to compute trends
38   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
39
40   !! * Substitutions
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
50      &                                       ptb, ptn, pta, kjpt )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE tra_adv_qck  ***
53      !!
54      !! ** Purpose :   Compute the now trend due to the advection of tracers
55      !!      and add it to the general trend of passive tracer equations.
56      !!
57      !! ** Method :   The advection is evaluated by a third order scheme
58      !!             For a positive velocity u :              u(i)>0
59      !!                                          |--FU--|--FC--|--FD--|------|
60      !!                                             i-1    i      i+1   i+2
61      !!
62      !!             For a negative velocity u :              u(i)<0
63      !!                                          |------|--FD--|--FC--|--FU--|
64      !!                                             i-1    i      i+1   i+2
65      !!             where  FU is the second upwind point
66      !!                    FD is the first douwning point
67      !!                    FC is the central point (or the first upwind point)
68      !!
69      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
70      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
71      !!
72      !!         dt = 2*rdtra and the scalar values are tb and sb
73      !!
74      !!       On the vertical, the simple centered scheme used ptn
75      !!
76      !!               The fluxes are bounded by the ULTIMATE limiter to
77      !!             guarantee the monotonicity of the solution and to
78      !!            prevent the appearance of spurious numerical oscillations
79      !!
80      !! ** Action : - update pta  with the now advective tracer trends
81      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
82      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
83      !!
84      !! ** Reference : Leonard (1979, 1991)
85      !!----------------------------------------------------------------------
86      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
87      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
88      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
89      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
90      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of tracer time-step
91      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
92      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
94      !!----------------------------------------------------------------------
95      !
96      IF( nn_timing == 1 )  CALL timing_start('tra_adv_qck')
97      !
98      IF( kt == kit000 )  THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
101         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
102         IF(lwp) WRITE(numout,*)
103      ENDIF
104      !
105      l_trd = .FALSE.
106      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
107      !
108      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
109      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
110      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
111
112      !        ! vertical fluxes are computed with the 2nd order centered scheme
113      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
114      !
115      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_qck')
116      !
117   END SUBROUTINE tra_adv_qck
118
119
120   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
121      &                                        ptb, ptn, pta, kjpt   )
122      !!----------------------------------------------------------------------
123      !!
124      !!----------------------------------------------------------------------
125      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
126      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
127      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
128      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
129      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-velocity components
130      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
131      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
132      !!
133      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
134      REAL(wp) ::   ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
135      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zfu, zfc, zfd
136      !----------------------------------------------------------------------
137      !
138      CALL wrk_alloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
139      !                                                          ! ===========
140      DO jn = 1, kjpt                                            ! tracer loop
141         !                                                       ! ===========
142         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
143         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
144         !
145!!gm why not using a SHIFT instruction...
146         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask
147            DO jj = 2, jpjm1
148               DO ji = fs_2, fs_jpim1   ! vector opt.
149                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
150                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
151               END DO
152            END DO
153         END DO
154         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
155         
156         !
157         ! Horizontal advective fluxes
158         ! ---------------------------
159         DO jk = 1, jpkm1                             
160            DO jj = 2, jpjm1
161               DO ji = fs_2, fs_jpim1   ! vector opt.         
162                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
163                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
164               END DO
165            END DO
166         END DO
167         !
168         DO jk = 1, jpkm1 
169            zdt =  p2dt(jk)
170            DO jj = 2, jpjm1
171               DO ji = fs_2, fs_jpim1   ! vector opt.   
172                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
173                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
174                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
175                  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
176                  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
177               END DO
178            END DO
179         END DO 
180         !--- Lateral boundary conditions
181         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
182         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
183
184         !--- QUICKEST scheme
185         CALL quickest( zfu, zfd, zfc, zwx )
186         !
187         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
188         DO jk = 1, jpkm1 
189            DO jj = 2, jpjm1
190               DO ji = fs_2, fs_jpim1   ! vector opt.               
191                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
192               END DO
193            END DO
194         END DO
195         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
196
197         !
198         ! Tracer flux on the x-direction
199         DO jk = 1, jpkm1 
200            !
201            DO jj = 2, jpjm1
202               DO ji = fs_2, fs_jpim1   ! vector opt.               
203                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
204                  !--- If the second ustream point is a land point
205                  !--- the flux is computed by the 1st order UPWIND scheme
206                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
207                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
208                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
209               END DO
210            END DO
211         END DO
212         !
213         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
214         !
215         ! Computation of the trend
216         DO jk = 1, jpkm1 
217            DO jj = 2, jpjm1
218               DO ji = fs_2, fs_jpim1   ! vector opt. 
219                  zbtr = 1. / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
220                  ! horizontal advective trends
221                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
222                  !--- add it to the general tracer trends
223                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
224               END DO
225            END DO
226         END DO
227         !                                 ! trend diagnostics
228         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
229         !
230      END DO
231      !
232      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
233      !
234   END SUBROUTINE tra_adv_qck_i
235
236
237   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
238      &                                        ptb, ptn, pta, kjpt )
239      !!----------------------------------------------------------------------
240      !!
241      !!----------------------------------------------------------------------
242      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
243      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
244      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
245      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
246      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
247      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
248      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
249      !!
250      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
251      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
252      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd
253      !----------------------------------------------------------------------
254      !
255      CALL wrk_alloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
256      !
257      !                                                          ! ===========
258      DO jn = 1, kjpt                                            ! tracer loop
259         !                                                       ! ===========
260         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
261         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
262         !                                                 
263         DO jk = 1, jpkm1                               
264            !                                             
265            !--- Computation of the ustream and downstream value of the tracer and the mask
266            DO jj = 2, jpjm1
267               DO ji = fs_2, fs_jpim1   ! vector opt.
268                  ! Upstream in the x-direction for the tracer
269                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
270                  ! Downstream in the x-direction for the tracer
271                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
272               END DO
273            END DO
274         END DO
275         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
276
277         
278         !
279         ! Horizontal advective fluxes
280         ! ---------------------------
281         !
282         DO jk = 1, jpkm1                             
283            DO jj = 2, jpjm1
284               DO ji = fs_2, fs_jpim1   ! vector opt.         
285                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
286                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
287               END DO
288            END DO
289         END DO
290         !
291         DO jk = 1, jpkm1 
292            zdt =  p2dt(jk)
293            DO jj = 2, jpjm1
294               DO ji = fs_2, fs_jpim1   ! vector opt.   
295                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
296                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
297                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
298                  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
299                  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
300               END DO
301            END DO
302         END DO
303
304         !--- Lateral boundary conditions
305         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
306         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
307
308         !--- QUICKEST scheme
309         CALL quickest( zfu, zfd, zfc, zwy )
310         !
311         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
312         DO jk = 1, jpkm1 
313            DO jj = 2, jpjm1
314               DO ji = fs_2, fs_jpim1   ! vector opt.               
315                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
316               END DO
317            END DO
318         END DO
319         !--- Lateral boundary conditions
320         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
321         !
322         ! Tracer flux on the x-direction
323         DO jk = 1, jpkm1 
324            !
325            DO jj = 2, jpjm1
326               DO ji = fs_2, fs_jpim1   ! vector opt.               
327                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
328                  !--- If the second ustream point is a land point
329                  !--- the flux is computed by the 1st order UPWIND scheme
330                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
331                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
332                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
333               END DO
334            END DO
335         END DO
336         !
337         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
338         !
339         ! Computation of the trend
340         DO jk = 1, jpkm1 
341            DO jj = 2, jpjm1
342               DO ji = fs_2, fs_jpim1   ! vector opt. 
343                  zbtr = 1. / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
344                  ! horizontal advective trends
345                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
346                  !--- add it to the general tracer trends
347                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
348               END DO
349            END DO
350         END DO
351         !                                 ! trend diagnostics
352         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
353         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
354         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
355           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
356           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
357         ENDIF
358         !
359      END DO
360      !
361      CALL wrk_dealloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
362      !
363   END SUBROUTINE tra_adv_qck_j
364
365
366   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
367     &                                    ptn, pta, kjpt )
368      !!----------------------------------------------------------------------
369      !!
370      !!----------------------------------------------------------------------
371      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
372      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
373      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
374      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
375      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
376      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
377      !
378      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
379      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz
380      !!----------------------------------------------------------------------
381      !
382      CALL wrk_alloc( jpi,jpj,jpk,   zwz )
383      !
384      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
385      zwz(:,:,jpk) = 0._wp
386      !
387      !                                                          ! ===========
388      DO jn = 1, kjpt                                            ! tracer loop
389         !                                                       ! ===========
390         !
391         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
392            DO jj = 2, jpjm1
393               DO ji = fs_2, fs_jpim1   ! vector opt.
394                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
395               END DO
396            END DO
397         END DO
398         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
399            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
400               DO jj = 1, jpj
401                  DO ji = 1, jpi
402                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
403                  END DO
404               END DO   
405            ELSE                                   ! no ocean cavities (only ocean surface)
406               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
407            ENDIF
408         ENDIF
409         !
410         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
411            DO jj = 2, jpjm1
412               DO ji = fs_2, fs_jpim1   ! vector opt.
413                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
414                     &                                / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
415               END DO
416            END DO
417         END DO
418         !                                 ! Send trends for diagnostic
419         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
420         !
421      END DO
422      !
423      CALL wrk_dealloc( jpi,jpj,jpk,   zwz )
424      !
425   END SUBROUTINE tra_adv_cen2_k
426
427
428   SUBROUTINE quickest( pfu, pfd, pfc, puc )
429      !!----------------------------------------------------------------------
430      !!
431      !! ** Purpose :  Computation of advective flux with Quickest scheme
432      !!
433      !! ** Method :   
434      !!----------------------------------------------------------------------
435      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
436      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
437      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
438      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
439      !!
440      INTEGER  ::  ji, jj, jk               ! dummy loop indices
441      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
442      REAL(wp) ::  zc, zcurv, zfho          !   -      -
443      !----------------------------------------------------------------------
444      !
445      IF( nn_timing == 1 )  CALL timing_start('quickest')
446      !
447      DO jk = 1, jpkm1
448         DO jj = 1, jpj
449            DO ji = 1, jpi
450               zc     = puc(ji,jj,jk)                         ! Courant number
451               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
452               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
453               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
454               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
455               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
456               !
457               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
458               zcoef2 = ABS( zcoef1 )
459               zcoef3 = ABS( zcurv )
460               IF( zcoef3 >= zcoef2 ) THEN
461                  zfho = pfc(ji,jj,jk) 
462               ELSE
463                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
464                  IF( zcoef1 >= 0. ) THEN
465                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
466                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
467                  ELSE
468                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
469                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
470                  ENDIF
471               ENDIF
472               puc(ji,jj,jk) = zfho
473            END DO
474         END DO
475      END DO
476      !
477      IF( nn_timing == 1 )  CALL timing_stop('quickest')
478      !
479   END SUBROUTINE quickest
480
481   !!======================================================================
482END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.