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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 5836

Last change on this file since 5836 was 5836, checked in by cetlod, 8 years ago

merge the simplification branch onto the trunk, see ticket #1612
  • Property svn:keywords set to Id
File size: 24.2 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 "domzgr_substitute.h90"
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
51      &                                       ptb, ptn, pta, kjpt )
52      !!----------------------------------------------------------------------
53      !!                  ***  ROUTINE tra_adv_qck  ***
54      !!
55      !! ** Purpose :   Compute the now trend due to the advection of tracers
56      !!      and add it to the general trend of passive tracer equations.
57      !!
58      !! ** Method :   The advection is evaluated by a third order scheme
59      !!             For a positive velocity u :              u(i)>0
60      !!                                          |--FU--|--FC--|--FD--|------|
61      !!                                             i-1    i      i+1   i+2
62      !!
63      !!             For a negative velocity u :              u(i)<0
64      !!                                          |------|--FD--|--FC--|--FU--|
65      !!                                             i-1    i      i+1   i+2
66      !!             where  FU is the second upwind point
67      !!                    FD is the first douwning point
68      !!                    FC is the central point (or the first upwind point)
69      !!
70      !!      Flux(i) = u(i) * { 0.5(FC+FD)  -0.5C(i)(FD-FC)  -((1-C(i))/6)(FU+FD-2FC) }
71      !!                with C(i)=|u(i)|dx(i)/dt (=Courant number)
72      !!
73      !!         dt = 2*rdtra and the scalar values are tb and sb
74      !!
75      !!       On the vertical, the simple centered scheme used ptn
76      !!
77      !!               The fluxes are bounded by the ULTIMATE limiter to
78      !!             guarantee the monotonicity of the solution and to
79      !!            prevent the appearance of spurious numerical oscillations
80      !!
81      !! ** Action : - update (pta) with the now advective tracer trends
82      !!             - save the trends
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      ! I. The 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      ! II. The 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) * fse3u(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) * fse3t(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 (contribution of upstream fluxes)
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) * fse3v(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) * fse3t(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 (contribution of upstream fluxes)
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      !                          ! surface & bottom values
385      IF( lk_vvl )   zwz(:,:, 1 ) = 0._wp             ! set to zero one for all
386                     zwz(:,:,jpk) = 0._wp             ! except at the surface in linear free surface
387      !
388      !                                                          ! ===========
389      DO jn = 1, kjpt                                            ! tracer loop
390         !                                                       ! ===========
391         !
392         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
393            DO jj = 2, jpjm1
394               DO ji = fs_2, fs_jpim1   ! vector opt.
395                  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)
396               END DO
397            END DO
398         END DO
399         IF(.NOT.lk_vvl ) THEN               !* top value   (only in linear free surf. as zwz is multiplied by wmask)
400            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
401               DO jj = 1, jpj
402                  DO ji = 1, jpi
403                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
404                  END DO
405               END DO   
406            ELSE                                   ! no ice-shelf cavities (only ocean surface)
407               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
408            ENDIF
409         ENDIF
410         !
411         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
412            DO jj = 2, jpjm1
413               DO ji = fs_2, fs_jpim1   ! vector opt.
414                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
415                     &                                / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
416               END DO
417            END DO
418         END DO
419         !                                 ! Save the vertical advective trends for diagnostic
420         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
421         !
422      END DO
423      !
424      CALL wrk_dealloc( jpi, jpj, jpk, zwz )
425      !
426   END SUBROUTINE tra_adv_cen2_k
427
428
429   SUBROUTINE quickest( pfu, pfd, pfc, puc )
430      !!----------------------------------------------------------------------
431      !!
432      !! ** Purpose :  Computation of advective flux with Quickest scheme
433      !!
434      !! ** Method :   
435      !!----------------------------------------------------------------------
436      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
437      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
438      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
439      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
440      !!
441      INTEGER  ::  ji, jj, jk               ! dummy loop indices
442      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
443      REAL(wp) ::  zc, zcurv, zfho          !   -      -
444      !----------------------------------------------------------------------
445      !
446      IF( nn_timing == 1 )  CALL timing_start('quickest')
447      !
448      DO jk = 1, jpkm1
449         DO jj = 1, jpj
450            DO ji = 1, jpi
451               zc     = puc(ji,jj,jk)                         ! Courant number
452               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
453               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
454               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
455               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
456               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
457               !
458               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
459               zcoef2 = ABS( zcoef1 )
460               zcoef3 = ABS( zcurv )
461               IF( zcoef3 >= zcoef2 ) THEN
462                  zfho = pfc(ji,jj,jk) 
463               ELSE
464                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
465                  IF( zcoef1 >= 0. ) THEN
466                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
467                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
468                  ELSE
469                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
470                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
471                  ENDIF
472               ENDIF
473               puc(ji,jj,jk) = zfho
474            END DO
475         END DO
476      END DO
477      !
478      IF( nn_timing == 1 )  CALL timing_stop('quickest')
479      !
480   END SUBROUTINE quickest
481
482   !!======================================================================
483END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.