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

source: branches/UKMO/dev_r8864_restart_date/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 9235

Last change on this file since 9235 was 9235, checked in by davestorkey, 6 years ago

UKMO/dev_r8864_restart_date : clear SVN keywords.

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   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
37
38   LOGICAL  ::   l_trd   ! flag to compute trends
39   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
40
41
42   !! * Substitutions
43#  include "vectopt_loop_substitute.h90"
44   !!----------------------------------------------------------------------
45   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
46   !! $Id$
47   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
48   !!----------------------------------------------------------------------
49CONTAINS
50
51   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
52      &                                       ptb, ptn, pta, kjpt )
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE tra_adv_qck  ***
55      !!
56      !! ** Purpose :   Compute the now trend due to the advection of tracers
57      !!      and add it to the general trend of passive tracer equations.
58      !!
59      !! ** Method :   The advection is evaluated by a third order scheme
60      !!             For a positive velocity u :              u(i)>0
61      !!                                          |--FU--|--FC--|--FD--|------|
62      !!                                             i-1    i      i+1   i+2
63      !!
64      !!             For a negative velocity u :              u(i)<0
65      !!                                          |------|--FD--|--FC--|--FU--|
66      !!                                             i-1    i      i+1   i+2
67      !!             where  FU is the second upwind point
68      !!                    FD is the first douwning point
69      !!                    FC is the central point (or the first upwind point)
70      !!
71      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
72      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
73      !!
74      !!         dt = 2*rdtra and the scalar values are tb and sb
75      !!
76      !!       On the vertical, the simple centered scheme used ptn
77      !!
78      !!               The fluxes are bounded by the ULTIMATE limiter to
79      !!             guarantee the monotonicity of the solution and to
80      !!            prevent the appearance of spurious numerical oscillations
81      !!
82      !! ** Action : - update pta  with the now advective tracer trends
83      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
84      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
85      !!
86      !! ** Reference : Leonard (1979, 1991)
87      !!----------------------------------------------------------------------
88      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
89      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
90      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
91      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
92      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
93      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
94      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
95      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
96      !!----------------------------------------------------------------------
97      !
98      IF( nn_timing == 1 )  CALL timing_start('tra_adv_qck')
99      !
100      IF( kt == kit000 )  THEN
101         IF(lwp) WRITE(numout,*)
102         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
103         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
104         IF(lwp) WRITE(numout,*)
105      ENDIF
106      !
107      l_trd = .FALSE.
108      l_ptr = .FALSE.
109      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
110      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
111      !
112      !
113      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
114      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
115      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
116
117      !        ! vertical fluxes are computed with the 2nd order centered scheme
118      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
119      !
120      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_qck')
121      !
122   END SUBROUTINE tra_adv_qck
123
124
125   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
126      &                                        ptb, ptn, pta, kjpt   )
127      !!----------------------------------------------------------------------
128      !!
129      !!----------------------------------------------------------------------
130      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
131      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
132      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
133      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
134      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun        ! i-velocity components
135      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
136      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
137      !!
138      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
139      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
140      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwx, zfu, zfc, zfd
141      !----------------------------------------------------------------------
142      !
143      CALL wrk_alloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
144      !                                                          ! ===========
145      DO jn = 1, kjpt                                            ! tracer loop
146         !                                                       ! ===========
147         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
148         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
149         !
150!!gm why not using a SHIFT instruction...
151         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask
152            DO jj = 2, jpjm1
153               DO ji = fs_2, fs_jpim1   ! vector opt.
154                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
155                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
156               END DO
157            END DO
158         END DO
159         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
160         
161         !
162         ! Horizontal advective fluxes
163         ! ---------------------------
164         DO jk = 1, jpkm1                             
165            DO jj = 2, jpjm1
166               DO ji = fs_2, fs_jpim1   ! vector opt.         
167                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
168                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
169               END DO
170            END DO
171         END DO
172         !
173         DO jk = 1, jpkm1 
174            DO jj = 2, jpjm1
175               DO ji = fs_2, fs_jpim1   ! vector opt.   
176                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
177                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
178                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
179                  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
180                  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
181               END DO
182            END DO
183         END DO 
184         !--- Lateral boundary conditions
185         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
186         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
187
188         !--- QUICKEST scheme
189         CALL quickest( zfu, zfd, zfc, zwx )
190         !
191         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
192         DO jk = 1, jpkm1 
193            DO jj = 2, jpjm1
194               DO ji = fs_2, fs_jpim1   ! vector opt.               
195                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
196               END DO
197            END DO
198         END DO
199         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
200
201         !
202         ! Tracer flux on the x-direction
203         DO jk = 1, jpkm1 
204            !
205            DO jj = 2, jpjm1
206               DO ji = fs_2, fs_jpim1   ! vector opt.               
207                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
208                  !--- If the second ustream point is a land point
209                  !--- the flux is computed by the 1st order UPWIND scheme
210                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
211                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
212                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
213               END DO
214            END DO
215         END DO
216         !
217         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
218         !
219         ! Computation of the trend
220         DO jk = 1, jpkm1 
221            DO jj = 2, jpjm1
222               DO ji = fs_2, fs_jpim1   ! vector opt. 
223                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
224                  ! horizontal advective trends
225                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
226                  !--- add it to the general tracer trends
227                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
228               END DO
229            END DO
230         END DO
231         !                                 ! trend diagnostics
232         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
233         !
234      END DO
235      !
236      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
237      !
238   END SUBROUTINE tra_adv_qck_i
239
240
241   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
242      &                                        ptb, ptn, pta, kjpt )
243      !!----------------------------------------------------------------------
244      !!
245      !!----------------------------------------------------------------------
246      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
247      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
248      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
249      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
250      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
251      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
252      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
253      !!
254      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
255      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
256      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd
257      !----------------------------------------------------------------------
258      !
259      CALL wrk_alloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
260      !
261      !                                                          ! ===========
262      DO jn = 1, kjpt                                            ! tracer loop
263         !                                                       ! ===========
264         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
265         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
266         !                                                 
267         DO jk = 1, jpkm1                               
268            !                                             
269            !--- Computation of the ustream and downstream value of the tracer and the mask
270            DO jj = 2, jpjm1
271               DO ji = fs_2, fs_jpim1   ! vector opt.
272                  ! Upstream in the x-direction for the tracer
273                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
274                  ! Downstream in the x-direction for the tracer
275                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
276               END DO
277            END DO
278         END DO
279         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
280
281         
282         !
283         ! Horizontal advective fluxes
284         ! ---------------------------
285         !
286         DO jk = 1, jpkm1                             
287            DO jj = 2, jpjm1
288               DO ji = fs_2, fs_jpim1   ! vector opt.         
289                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
290                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
291               END DO
292            END DO
293         END DO
294         !
295         DO jk = 1, jpkm1 
296            DO jj = 2, jpjm1
297               DO ji = fs_2, fs_jpim1   ! vector opt.   
298                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
299                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
300                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
301                  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
302                  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
303               END DO
304            END DO
305         END DO
306
307         !--- Lateral boundary conditions
308         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
309         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
310
311         !--- QUICKEST scheme
312         CALL quickest( zfu, zfd, zfc, zwy )
313         !
314         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
315         DO jk = 1, jpkm1 
316            DO jj = 2, jpjm1
317               DO ji = fs_2, fs_jpim1   ! vector opt.               
318                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
319               END DO
320            END DO
321         END DO
322         !--- Lateral boundary conditions
323         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
324         !
325         ! Tracer flux on the x-direction
326         DO jk = 1, jpkm1 
327            !
328            DO jj = 2, jpjm1
329               DO ji = fs_2, fs_jpim1   ! vector opt.               
330                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
331                  !--- If the second ustream point is a land point
332                  !--- the flux is computed by the 1st order UPWIND scheme
333                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
334                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
335                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
336               END DO
337            END DO
338         END DO
339         !
340         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
341         !
342         ! Computation of the trend
343         DO jk = 1, jpkm1 
344            DO jj = 2, jpjm1
345               DO ji = fs_2, fs_jpim1   ! vector opt. 
346                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
347                  ! horizontal advective trends
348                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
349                  !--- add it to the general tracer trends
350                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
351               END DO
352            END DO
353         END DO
354         !                                 ! trend diagnostics
355         IF( l_trd )                     CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
356         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
357         IF( l_ptr )                     CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
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                     &                                * r1_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.