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 NEMO/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/traadv_qck.F90 @ 11993

Last change on this file since 11993 was 11993, checked in by cetlod, 4 years ago

trunk : undo bad commit. Oups

  • Property svn:keywords set to Id
File size: 23.2 KB
RevLine 
[1231]1MODULE traadv_qck
2   !!==============================================================================
3   !!                       ***  MODULE  traadv_qck  ***
[2528]4   !! Ocean tracers:  horizontal & vertical advective trend
[1231]5   !!==============================================================================
[1559]6   !! History :  3.0  !  2008-07  (G. Reffray)  Original code
[2528]7   !!            3.3  !  2010-05  (C.Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
[1231]8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
[2528]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
[1559]15   !!   tra_adv_cen2_k : 2nd centered scheme for the vertical advection
[1231]16   !!----------------------------------------------------------------------
17   USE oce             ! ocean dynamics and active tracers
18   USE dom_oce         ! ocean space and time domain
[4990]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   !
[9124]24   USE in_out_manager  ! I/O manager
[1231]25   USE lib_mpp         ! distribued memory computing
26   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
[3625]27   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
[1231]28
29   IMPLICIT NONE
30   PRIVATE
31
[1559]32   PUBLIC   tra_adv_qck   ! routine called by step.F90
[1231]33
[2528]34   REAL(wp) :: r1_6 = 1./ 6.   ! 1/6 ratio
[1559]35
[7646]36   LOGICAL  ::   l_trd   ! flag to compute trends
37   LOGICAL  ::   l_ptr   ! flag to compute poleward transport
38
39
[1231]40   !! * Substitutions
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
[9598]43   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1231]44   !! $Id$
[10068]45   !! Software governed by the CeCILL license (see ./LICENSE)
[1231]46   !!----------------------------------------------------------------------
47CONTAINS
48
[3294]49   SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,      &
[2528]50      &                                       ptb, ptn, pta, kjpt )
[1231]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
[1559]58      !!             For a positive velocity u :              u(i)>0
59      !!                                          |--FU--|--FC--|--FD--|------|
60      !!                                             i-1    i      i+1   i+2
[1231]61      !!
[1559]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)
[1231]68      !!
[1559]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)
[1231]71      !!
72      !!         dt = 2*rdtra and the scalar values are tb and sb
73      !!
[2528]74      !!       On the vertical, the simple centered scheme used ptn
[1231]75      !!
[1559]76      !!               The fluxes are bounded by the ULTIMATE limiter to
77      !!             guarantee the monotonicity of the solution and to
[1231]78      !!            prevent the appearance of spurious numerical oscillations
79      !!
[6140]80      !! ** Action : - update pta  with the now advective tracer trends
81      !!             - send trends to trdtra module for further diagnostcs (l_trdtra=T)
[11993]82      !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
[1231]83      !!
84      !! ** Reference : Leonard (1979, 1991)
85      !!----------------------------------------------------------------------
[2528]86      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
[3294]87      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
[2528]88      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
89      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
[6140]90      REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step
[2528]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
[1231]94      !!----------------------------------------------------------------------
[3294]95      !
96      IF( kt == kit000 )  THEN
[1231]97         IF(lwp) WRITE(numout,*)
[2528]98         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
[1231]99         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
100         IF(lwp) WRITE(numout,*)
101      ENDIF
[5836]102      !
[4990]103      l_trd = .FALSE.
[7646]104      l_ptr = .FALSE.
[11993]105      IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )      l_trd = .TRUE.
106      IF(   cdtype == 'TRA' .AND. ln_diaptr )                                               l_ptr = .TRUE. 
[4499]107      !
[7646]108      !
[6140]109      !        ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
[2528]110      CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt ) 
111      CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt ) 
[1231]112
[6140]113      !        ! vertical fluxes are computed with the 2nd order centered scheme
[2528]114      CALL tra_adv_cen2_k( kt, cdtype, pwn,         ptn, pta, kjpt )
[1231]115      !
116   END SUBROUTINE tra_adv_qck
117
118
[2528]119   SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun,                  &
120      &                                        ptb, ptn, pta, kjpt   )
[1231]121      !!----------------------------------------------------------------------
122      !!
123      !!----------------------------------------------------------------------
[2715]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
[6140]127      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
[2715]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
[2528]131      !!
[5836]132      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
[6140]133      REAL(wp) ::   ztra, zbtr, zdir, zdx, zmsk   ! local scalars
[9019]134      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zfu, zfc, zfd
[1231]135      !----------------------------------------------------------------------
[2715]136      !
[2528]137      !                                                          ! ===========
138      DO jn = 1, kjpt                                            ! tracer loop
139         !                                                       ! ===========
[5836]140         zfu(:,:,:) = 0._wp     ;   zfc(:,:,:) = 0._wp 
141         zfd(:,:,:) = 0._wp     ;   zwx(:,:,:) = 0._wp   
142         !
143!!gm why not using a SHIFT instruction...
144         DO jk = 1, jpkm1     !--- Computation of the ustream and downstream value of the tracer and the mask
[2528]145            DO jj = 2, jpjm1
146               DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]147                  zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)        ! Upstream   in the x-direction for the tracer
148                  zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)        ! Downstream in the x-direction for the tracer
[2528]149               END DO
[1559]150            END DO
151         END DO
[10425]152         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
[2528]153         
[1231]154         !
155         ! Horizontal advective fluxes
156         ! ---------------------------
[2528]157         DO jk = 1, jpkm1                             
158            DO jj = 2, jpjm1
159               DO ji = fs_2, fs_jpim1   ! vector opt.         
160                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
161                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk)  ! FU in the x-direction for T
162               END DO
163            END DO
[1559]164         END DO
[1231]165         !
[2528]166         DO jk = 1, jpkm1 
167            DO jj = 2, jpjm1
168               DO ji = fs_2, fs_jpim1   ! vector opt.   
169                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
[6140]170                  zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
171                  zwx(ji,jj,jk)  = ABS( pun(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
[2528]172                  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
173                  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
174               END DO
175            END DO
176         END DO 
177         !--- Lateral boundary conditions
[10425]178         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1.,  zwx(:,:,:), 'T', 1. )
[2528]179
[1231]180         !--- QUICKEST scheme
[2528]181         CALL quickest( zfu, zfd, zfc, zwx )
[1231]182         !
[2528]183         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
184         DO jk = 1, jpkm1 
185            DO jj = 2, jpjm1
186               DO ji = fs_2, fs_jpim1   ! vector opt.               
187                  zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
[2715]188               END DO
[1231]189            END DO
190         END DO
[10425]191         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )      ! Lateral boundary conditions
[2528]192
[1231]193         !
[2528]194         ! Tracer flux on the x-direction
195         DO jk = 1, jpkm1 
196            !
[1231]197            DO jj = 2, jpjm1
[2528]198               DO ji = fs_2, fs_jpim1   ! vector opt.               
199                  zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
200                  !--- If the second ustream point is a land point
201                  !--- the flux is computed by the 1st order UPWIND scheme
202                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
203                  zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
204                  zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
[1231]205               END DO
206            END DO
[3300]207         END DO
208         !
[10425]209         CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
[3300]210         !
211         ! Computation of the trend
212         DO jk = 1, jpkm1 
[2528]213            DO jj = 2, jpjm1
214               DO ji = fs_2, fs_jpim1   ! vector opt. 
[6140]215                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[2528]216                  ! horizontal advective trends
217                  ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
218                  !--- add it to the general tracer trends
219                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
220               END DO
221            END DO
[1231]222         END DO
[6140]223         !                                 ! trend diagnostics
[9019]224         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
[2528]225         !
226      END DO
227      !
[1559]228   END SUBROUTINE tra_adv_qck_i
[1231]229
230
[2528]231   SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn,                &
232      &                                        ptb, ptn, pta, kjpt )
[1231]233      !!----------------------------------------------------------------------
234      !!
235      !!----------------------------------------------------------------------
[2715]236      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
237      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
238      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
[6140]239      REAL(wp)                             , INTENT(in   ) ::   p2dt       ! tracer time-step
[2715]240      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pvn        ! j-velocity components
241      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn   ! before and now tracer fields
242      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
[1559]243      !!
[9019]244      INTEGER  :: ji, jj, jk, jn                ! dummy loop indices
[6140]245      REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk   ! local scalars
[9019]246      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, zfu, zfc, zfd   ! 3D workspace
[1231]247      !----------------------------------------------------------------------
[2715]248      !
[2528]249      !                                                          ! ===========
250      DO jn = 1, kjpt                                            ! tracer loop
251         !                                                       ! ===========
252         zfu(:,:,:) = 0.0     ;   zfc(:,:,:) = 0.0 
253         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0     
254         !                                                 
255         DO jk = 1, jpkm1                               
256            !                                             
257            !--- Computation of the ustream and downstream value of the tracer and the mask
258            DO jj = 2, jpjm1
259               DO ji = fs_2, fs_jpim1   ! vector opt.
260                  ! Upstream in the x-direction for the tracer
261                  zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
262                  ! Downstream in the x-direction for the tracer
263                  zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
264               END DO
[1559]265            END DO
266         END DO
[10425]267         CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. )   ! Lateral boundary conditions
[2528]268
269         
[1231]270         !
271         ! Horizontal advective fluxes
272         ! ---------------------------
273         !
[2528]274         DO jk = 1, jpkm1                             
275            DO jj = 2, jpjm1
276               DO ji = fs_2, fs_jpim1   ! vector opt.         
277                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
278                  zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk)  ! FU in the x-direction for T
279               END DO
[1559]280            END DO
281         END DO
[1231]282         !
[2528]283         DO jk = 1, jpkm1 
284            DO jj = 2, jpjm1
285               DO ji = fs_2, fs_jpim1   ! vector opt.   
286                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
[6140]287                  zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
288                  zwy(ji,jj,jk)  = ABS( pvn(ji,jj,jk) ) * p2dt / zdx    ! (0<zc_cfl<1 : Courant number on x-direction)
[2528]289                  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
290                  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
291               END DO
292            END DO
293         END DO
294
295         !--- Lateral boundary conditions
[10425]296         CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. )
[2528]297
[1231]298         !--- QUICKEST scheme
[2528]299         CALL quickest( zfu, zfd, zfc, zwy )
[1231]300         !
[2528]301         ! Mask at the T-points in the x-direction (mask=0 or mask=1)
302         DO jk = 1, jpkm1 
303            DO jj = 2, jpjm1
304               DO ji = fs_2, fs_jpim1   ! vector opt.               
305                  zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
306               END DO
[1231]307            END DO
308         END DO
[10425]309         CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. )    !--- Lateral boundary conditions
[2528]310         !
311         ! Tracer flux on the x-direction
312         DO jk = 1, jpkm1 
313            !
[1231]314            DO jj = 2, jpjm1
[2528]315               DO ji = fs_2, fs_jpim1   ! vector opt.               
316                  zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )   ! if pun > 0 : zdir = 1 otherwise zdir = 0
317                  !--- If the second ustream point is a land point
318                  !--- the flux is computed by the 1st order UPWIND scheme
319                  zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
320                  zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
321                  zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
[1231]322               END DO
323            END DO
[3300]324         END DO
325         !
[10425]326         CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
[3300]327         !
328         ! Computation of the trend
329         DO jk = 1, jpkm1 
[2528]330            DO jj = 2, jpjm1
331               DO ji = fs_2, fs_jpim1   ! vector opt. 
[6140]332                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[2528]333                  ! horizontal advective trends
334                  ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
335                  !--- add it to the general tracer trends
336                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[1231]337               END DO
338            END DO
[2528]339         END DO
[6140]340         !                                 ! trend diagnostics
[9019]341         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
[2528]342         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[9019]343         IF( l_ptr )   CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
[2528]344         !
345      END DO
346      !
[1559]347   END SUBROUTINE tra_adv_qck_j
[1231]348
349
[2528]350   SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn,           &
351     &                                    ptn, pta, kjpt )
[1231]352      !!----------------------------------------------------------------------
353      !!
354      !!----------------------------------------------------------------------
[2715]355      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
356      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
357      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
358      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pwn      ! vertical velocity
359      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptn      ! before and now tracer fields
360      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
361      !
[2528]362      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
[9019]363      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz   ! 3D workspace
[1559]364      !!----------------------------------------------------------------------
[4990]365      !
[6140]366      zwz(:,:, 1 ) = 0._wp       ! surface & bottom values set to zero for all tracers
367      zwz(:,:,jpk) = 0._wp
[5836]368      !
[2528]369      !                                                          ! ===========
370      DO jn = 1, kjpt                                            ! tracer loop
371         !                                                       ! ===========
372         !
[5836]373         DO jk = 2, jpkm1                    !* Interior point   (w-masked 2nd order centered flux)
[2528]374            DO jj = 2, jpjm1
375               DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]376                  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)
[2528]377               END DO
[1231]378            END DO
379         END DO
[6140]380         IF( ln_linssh ) THEN                !* top value   (only in linear free surf. as zwz is multiplied by wmask)
[5836]381            IF( ln_isfcav ) THEN                  ! ice-shelf cavities (top of the ocean)
382               DO jj = 1, jpj
383                  DO ji = 1, jpi
384                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn)   ! linear free surface
385                  END DO
386               END DO   
[6140]387            ELSE                                   ! no ocean cavities (only ocean surface)
[5836]388               zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
389            ENDIF
390         ENDIF
[2528]391         !
392         DO jk = 1, jpkm1          !==  Tracer flux divergence added to the general trend  ==!
393            DO jj = 2, jpjm1
394               DO ji = fs_2, fs_jpim1   ! vector opt.
[5836]395                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )   &
[6140]396                     &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[2528]397               END DO
[1231]398            END DO
399         END DO
[6140]400         !                                 ! Send trends for diagnostic
[4990]401         IF( l_trd )  CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
[2528]402         !
[1231]403      END DO
404      !
[1559]405   END SUBROUTINE tra_adv_cen2_k
[1231]406
407
[2528]408   SUBROUTINE quickest( pfu, pfd, pfc, puc )
[1231]409      !!----------------------------------------------------------------------
410      !!
[2528]411      !! ** Purpose :  Computation of advective flux with Quickest scheme
412      !!
413      !! ** Method :   
[1231]414      !!----------------------------------------------------------------------
[2528]415      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfu   ! second upwind point
416      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfd   ! first douwning point
417      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pfc   ! the central point (or the first upwind point)
418      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   puc   ! input as Courant number ; output as flux
419      !!
420      INTEGER  ::  ji, jj, jk               ! dummy loop indices
421      REAL(wp) ::  zcoef1, zcoef2, zcoef3   ! local scalars         
422      REAL(wp) ::  zc, zcurv, zfho          !   -      -
423      !----------------------------------------------------------------------
[3294]424      !
[2528]425      DO jk = 1, jpkm1
426         DO jj = 1, jpj
427            DO ji = 1, jpi
428               zc     = puc(ji,jj,jk)                         ! Courant number
429               zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
430               zcoef1 = 0.5 *      ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
431               zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
432               zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
433               zfho   = zcoef1 - zcoef2 - zcoef3              !  phi_f QUICKEST
434               !
435               zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
436               zcoef2 = ABS( zcoef1 )
437               zcoef3 = ABS( zcurv )
438               IF( zcoef3 >= zcoef2 ) THEN
439                  zfho = pfc(ji,jj,jk) 
440               ELSE
441                  zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) )    ! phi_REF
442                  IF( zcoef1 >= 0. ) THEN
443                     zfho = MAX( pfc(ji,jj,jk), zfho ) 
444                     zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) ) 
445                  ELSE
446                     zfho = MIN( pfc(ji,jj,jk), zfho ) 
447                     zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) ) 
448                  ENDIF
449               ENDIF
450               puc(ji,jj,jk) = zfho
451            END DO
452         END DO
453      END DO
[1231]454      !
455   END SUBROUTINE quickest
456
457   !!======================================================================
458END MODULE traadv_qck
Note: See TracBrowser for help on using the repository browser.