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

source: branches/2015/dev_CMCC_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 6051

Last change on this file since 6051 was 6051, checked in by lovato, 8 years ago

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