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

source: branches/UKMO/r5936_hadgem3_cplfld/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 7138

Last change on this file since 7138 was 7138, checked in by jcastill, 7 years ago

Remove svn keywords

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 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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt            ! vertical profile of 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), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of 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, zdt, 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            zdt =  p2dt(jk)
169            DO jj = 2, jpjm1
170               DO ji = fs_2, fs_jpim1   ! vector opt.   
171                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
172                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
173                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
174                  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
175                  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
176               END DO
177            END DO
178         END DO 
179         !--- Lateral boundary conditions
180         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
181         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
182
183         !--- QUICKEST scheme
184         CALL quickest( zfu, zfd, zfc, zwx )
185         !
186         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
187         DO jk = 1, jpkm1 
188            DO jj = 2, jpjm1
189               DO ji = fs_2, fs_jpim1   ! vector opt.               
190                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
191               END DO
192            END DO
193         END DO
194         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
195
196         !
197         ! Tracer flux on the x-direction
198         DO jk = 1, jpkm1 
199            !
200            DO jj = 2, jpjm1
201               DO ji = fs_2, fs_jpim1   ! vector opt.               
202                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
203                  !--- If the second ustream point is a land point
204                  !--- the flux is computed by the 1st order UPWIND scheme
205                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
206                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
207                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
208               END DO
209            END DO
210         END DO
211         !
212         CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
213         !
214         ! Computation of the trend
215         DO jk = 1, jpkm1 
216            DO jj = 2, jpjm1
217               DO ji = fs_2, fs_jpim1   ! vector opt. 
218                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
219                  ! horizontal advective trends
220                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
221                  !--- add it to the general tracer trends
222                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
223               END DO
224            END DO
225         END DO
226         !                                 ! trend diagnostics (contribution of upstream fluxes)
227         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
228         !
229      END DO
230      !
231      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
232      !
233   END SUBROUTINE tra_adv_qck_i
234
235
236   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
237      &                                        ptb, ptn, pta, kjpt )
238      !!----------------------------------------------------------------------
239      !!
240      !!----------------------------------------------------------------------
241      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
242      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
243      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
244      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt       ! vertical profile of tracer time-step
245      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
246      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
247      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
248      !!
249      INTEGER  :: ji, jj, jk, jn   ! dummy loop indices
250      REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk   ! local scalars
251      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd
252      !----------------------------------------------------------------------
253      !
254      CALL wrk_alloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
255      !
256      !                                                          ! ===========
257      DO jn = 1, kjpt                                            ! tracer loop
258         !                                                       ! ===========
259         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
260         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
261         !                                                 
262         DO jk = 1, jpkm1                               
263            !                                             
264            !--- Computation of the ustream and downstream value of the tracer and the mask
265            DO jj = 2, jpjm1
266               DO ji = fs_2, fs_jpim1   ! vector opt.
267                  ! Upstream in the x-direction for the tracer
268                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
269                  ! Downstream in the x-direction for the tracer
270                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
271               END DO
272            END DO
273         END DO
274         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
275
276         
277         !
278         ! Horizontal advective fluxes
279         ! ---------------------------
280         !
281         DO jk = 1, jpkm1                             
282            DO jj = 2, jpjm1
283               DO ji = fs_2, fs_jpim1   ! vector opt.         
284                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
285                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
286               END DO
287            END DO
288         END DO
289         !
290         DO jk = 1, jpkm1 
291            zdt =  p2dt(jk)
292            DO jj = 2, jpjm1
293               DO ji = fs_2, fs_jpim1   ! vector opt.   
294                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
295                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
296                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * zdt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
297                  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
298                  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
299               END DO
300            END DO
301         END DO
302
303         !--- Lateral boundary conditions
304         CALL lbc_lnk( zfu(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
305         CALL lbc_lnk( zfc(:,:,:), 'T', 1. )      ;     CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
306
307         !--- QUICKEST scheme
308         CALL quickest( zfu, zfd, zfc, zwy )
309         !
310         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
311         DO jk = 1, jpkm1 
312            DO jj = 2, jpjm1
313               DO ji = fs_2, fs_jpim1   ! vector opt.               
314                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
315               END DO
316            END DO
317         END DO
318         !--- Lateral boundary conditions
319         CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) 
320         !
321         ! Tracer flux on the x-direction
322         DO jk = 1, jpkm1 
323            !
324            DO jj = 2, jpjm1
325               DO ji = fs_2, fs_jpim1   ! vector opt.               
326                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
327                  !--- If the second ustream point is a land point
328                  !--- the flux is computed by the 1st order UPWIND scheme
329                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
330                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
331                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
332               END DO
333            END DO
334         END DO
335         !
336         CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
337         !
338         ! Computation of the trend
339         DO jk = 1, jpkm1 
340            DO jj = 2, jpjm1
341               DO ji = fs_2, fs_jpim1   ! vector opt. 
342                  zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
343                  ! horizontal advective trends
344                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
345                  !--- add it to the general tracer trends
346                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
347               END DO
348            END DO
349         END DO
350         !                                 ! trend diagnostics (contribution of upstream fluxes)
351         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
352         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
353         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
354           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
355           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
356         ENDIF
357         !
358      END DO
359      !
360      CALL wrk_dealloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
361      !
362   END SUBROUTINE tra_adv_qck_j
363
364
365   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
366     &                                    ptn, pta, kjpt )
367      !!----------------------------------------------------------------------
368      !!
369      !!----------------------------------------------------------------------
370      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
371      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
372      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
373      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
374      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
375      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
376      !
377      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
378      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz
379      !!----------------------------------------------------------------------
380      !
381      CALL wrk_alloc( jpi,jpj,jpk,   zwz )
382      !
383      !                          ! surface & bottom values
384      IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp             ! set to zero one for all
385                     zwz(:,:,jpk) = 0._wp             ! except at the surface in linear free surface
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(.NOT.lk_vvl ) 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 ice-shelf 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) * fse3t(ji,jj,jk) )
415               END DO
416            END DO
417         END DO
418         !                                 ! Save the vertical advective 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.