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.
Changeset 791 – NEMO

Changeset 791


Ignore:
Timestamp:
2008-01-12T21:33:34+01:00 (16 years ago)
Author:
gm
Message:

dev_001_GM - TRA/traadv : switch from velocity to transport + optimised traadv_eiv2 introduced - compilation OK

Location:
branches/dev_001_GM/NEMO/OPA_SRC/TRA
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv.F90

    r786 r791  
    44   !! Ocean active tracers:  advection trend  
    55   !!============================================================================== 
    6    !! History :  9.0  !  05-11  (G. Madec)  Original code 
     6   !! History :  1.0  !  2005-11  (G. Madec)  Original code 
     7   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport 
    78   !!---------------------------------------------------------------------- 
    89 
     
    2324   USE ldftra_oce      ! lateral diffusion coefficient on tracers 
    2425   USE in_out_manager  ! I/O manager 
    25 !  USE prtctl          ! Print control 
    2626 
    2727   IMPLICIT NONE 
     
    5959      !! ** Method  : - Update (ua,va) with the advection term following nadv 
    6060      !!---------------------------------------------------------------------- 
    61 #if ( defined key_trabbl_adv || defined key_traldf_eiv ) 
    62       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zun, zvn, zwn   ! effective velocity 
    63 #else 
    64       USE oce, ONLY :                       zun => un       ! the effective velocity is the 
    65       USE oce, ONLY :                       zvn => vn       ! Eulerian velocity 
    66       USE oce, ONLY :                       zwn => wn       !  
    67 #endif 
     61      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    6862      !! 
    69       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
     63      INTEGER ::   jk   ! dummy loop index 
     64      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zun, zvn, zwn   ! effective transports 
    7065      !!---------------------------------------------------------------------- 
    7166 
    7267      IF( kt == nit000 )   CALL tra_adv_ctl          ! initialisation & control of options 
    7368 
     69 
     70      !                                              ! effective transport 
     71      DO jk = 1, jpkm1  
    7472#if defined key_trabbl_adv 
    75       zun(:,:,:) = un(:,:,:) - u_bbl(:,:,:)          ! add the bbl velocity 
    76       zvn(:,:,:) = vn(:,:,:) - v_bbl(:,:,:) 
    77       zwn(:,:,:) = wn(:,:,:) + w_bbl(:,:,:) 
     73         !                                                ! eulerian + bbl transport 
     74         zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * ( un(:,:,jk) - u_bbl(:,:,jk) ) 
     75         zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * ( vn(:,:,jk) - v_bbl(:,:,jk) ) 
     76         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * ( wn(:,:,jk) + w_bbl(:,:,jk) ) 
     77#else 
     78         !                                                ! eulerian transport only 
     79         zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) *   un(:,:,jk)  
     80         zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) *   vn(:,:,jk) 
     81         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      *   wn(:,:,jk) 
    7882#endif 
    79       IF( lk_traldf_eiv ) THEN                       ! commpute and add the eiv velocity 
    80          IF( .NOT. lk_trabbl_adv ) THEN  
    81             zun(:,:,:) = un(:,:,:) 
    82             zvn(:,:,:) = vn(:,:,:) 
    83             zwn(:,:,:) = wn(:,:,:) 
    84          ENDIF 
    85          CALL tra_adv_eiv( kt, zun, zvn, zwn )  
    86       ENDIF 
     83      END DO 
     84      zwn(:,:,jpk) = 0.e0                                 ! no transport trough the bottom 
     85 
     86      !                                                   ! add the eiv transport (if necessary) 
     87      IF( lk_traldf_eiv )   CALL tra_adv_eiv( kt, zun, zvn, zwn )  
     88 
    8789 
    8890      SELECT CASE ( nadv )                           ! compute advection trend and add it to general trend 
    89       CASE ( 0 )   ;   CALL tra_adv_cen2    ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! 2nd order centered 
    90                        CALL tra_adv_cen2    ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! 2nd order centered 
    91 !     CASE ( 1 )   ;   CALL tra_adv_cen2_jki( kt, zun, zvn, zwn )    ! 2nd order centered scheme 
     91      ! 
     92      CASE ( 1 )   ;   CALL tra_adv_cen2    ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! 2nd order centered 
     93                       CALL tra_adv_cen2    ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S 
     94      ! 
    9295      CASE ( 2 )   ;   CALL tra_adv_tvd     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! TVD      scheme 
    93                        CALL tra_adv_tvd     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! TVD      scheme 
     96                       CALL tra_adv_tvd     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S 
     97      ! 
    9498      CASE ( 3 )   ;   CALL tra_adv_muscl   ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb    , ta )    ! MUSCL    scheme 
    95                        CALL tra_adv_muscl   ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb    , sa )    ! MUSCL    scheme 
     99                       CALL tra_adv_muscl   ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb    , sa )    ! on T & S 
     100      ! 
    96101      CASE ( 4 )   ;   CALL tra_adv_muscl2  ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! MUSCL2   scheme 
    97                        CALL tra_adv_muscl2  ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! MUSCL2   scheme 
     102                       CALL tra_adv_muscl2  ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S 
     103      ! 
    98104      CASE ( 5 )   ;   CALL tra_adv_ubs     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! UBS      scheme 
    99                        CALL tra_adv_ubs     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! UBS      scheme 
     105                       CALL tra_adv_ubs     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S 
     106      ! 
    100107      CASE ( 6 )   ;   CALL tra_adv_qck     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! QUICKEST scheme 
    101                        CALL tra_adv_qck     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! QUICKEST scheme 
     108                       CALL tra_adv_qck     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! on T & S 
    102109      ! 
    103       CASE (-1 )                                                     ! esopa: test all possibility with control print 
     110      CASE (-1 )                                     ! NEMO debug : test all the schemes at once 
    104111                       CALL tra_adv_cen2    ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! 2nd order centered 
    105                        CALL tra_adv_cen2    ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! 2nd order centered 
     112                       CALL tra_adv_cen2    ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )  
    106113                       CALL tra_adv_tvd     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! TVD      scheme 
    107                        CALL tra_adv_tvd     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! TVD      scheme 
     114                       CALL tra_adv_tvd     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )   
    108115                       CALL tra_adv_muscl   ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb    , ta )    ! MUSCL    scheme 
    109                        CALL tra_adv_muscl   ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb    , sa )    ! MUSCL    scheme 
     116                       CALL tra_adv_muscl   ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb    , sa )  
    110117                       CALL tra_adv_muscl2  ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! MUSCL2   scheme 
    111                        CALL tra_adv_muscl2  ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! MUSCL2   scheme 
     118                       CALL tra_adv_muscl2  ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) 
    112119                       CALL tra_adv_ubs     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! UBS      scheme 
    113                        CALL tra_adv_ubs     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! UBS      scheme 
     120                       CALL tra_adv_ubs     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) 
    114121                       CALL tra_adv_qck     ( kt, 'TRA', jp_tem, zun, zvn, zwn, tb, tn, ta )    ! QUICKEST scheme 
    115                        CALL tra_adv_qck     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa )    ! QUICKEST scheme 
     122                       CALL tra_adv_qck     ( kt, 'TRA', jp_sal, zun, zvn, zwn, sb, sn, sa ) 
    116123      END SELECT 
    117124      ! 
     
    127134      !!---------------------------------------------------------------------- 
    128135      INTEGER ::   ioptio 
    129  
     136      !! 
    130137      NAMELIST/nam_traadv/ ln_traadv_cen2 , ln_traadv_tvd,    & 
    131138         &                 ln_traadv_muscl, ln_traadv_muscl2, & 
     
    163170         &                CALL ctl_stop( 'cross-land advection only with 2nd order advection scheme' ) 
    164171 
    165       !                              ! Set nadv 
    166       IF( ln_traadv_cen2   )   nadv =  0 
    167 #if defined key_mpp_omp 
     172      !                               ! Set nadv 
    168173      IF( ln_traadv_cen2   )   nadv =  1 
    169 #endif 
    170174      IF( ln_traadv_tvd    )   nadv =  2 
    171175      IF( ln_traadv_muscl  )   nadv =  3 
     
    175179      IF( lk_esopa         )   nadv = -1 
    176180 
    177       IF(lwp) THEN                   ! Print the choice 
     181      IF(lwp) THEN                    ! Print the choice 
    178182         WRITE(numout,*) 
    179          IF( nadv ==  0 )   WRITE(numout,*) '         2nd order scheme is used' 
    180          IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is usedi, k-j-i case' 
     183         IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used' 
    181184         IF( nadv ==  2 )   WRITE(numout,*) '         TVD       scheme is used' 
    182185         IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     scheme is used' 
     
    185188         IF( nadv ==  6 )   WRITE(numout,*) '         PPM       scheme is used' 
    186189         IF( nadv ==  7 )   WRITE(numout,*) '         QUICKEST  scheme is used' 
    187          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme' 
     190         IF( nadv == -1 )   WRITE(numout,*) '         NEMO debug: test all advection scheme at once' 
    188191      ENDIF 
    189192      ! 
  • branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90

    r786 r791  
    44   !! Ocean active tracers:  horizontal & vertical advective trend 
    55   !!============================================================================== 
    6    !! History :  8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv  
    7    !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
    8    !!   NEMO     1.0  !  05-11  (V. Garnier) Surface pressure gradient organization 
    9    !!             -   !  05-11  (V. Garnier) Surface pressure gradient organization 
    10    !!             -   !  06-04  (R. Benshila, G. Madec) Step reorganization 
    11    !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC 
     6   !! History :  8.2  !  2001-08  (G. Madec, E. Durand)  trahad+trazad = traadv  
     7   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module 
     8   !!   NEMO     1.0  !  2005-11  (V. Garnier) Surface pressure gradient organization 
     9   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization 
     10   !!             -   !  2006-04  (R. Benshila, G. Madec) Step reorganization 
     11   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport 
    1212   !!---------------------------------------------------------------------- 
    1313 
     
    3434   PUBLIC   tra_adv_cen2   ! routine called by step.F90 
    3535 
    36    REAL(wp), DIMENSION(jpi,jpj) ::   btr2   ! inverse of T-point surface [1/(e1t*e2t)] 
    37  
    3836   !! * Substitutions 
    3937#  include "domzgr_substitute.h90" 
     
    4139   !!---------------------------------------------------------------------- 
    4240   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)  
    43    !! $Id:$  
     41   !! $Id$  
    4442   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4543   !!---------------------------------------------------------------------- 
     
    6260      !!         Part 0 : compute the upstream / centered flag 
    6361      !!                  (3D array, zind, defined at T-point (0<zind<1)) 
    64       !!         Part I : horizontal advection 
    65       !!       * centered flux: 
    66       !!               zcenu = e2u*e3u  un  mi(ptn) 
    67       !!               zcenv = e1v*e3v  vn  mj(ptn) 
    68       !!       * upstream flux: 
    69       !!               zupsu = e2u*e3u  un  (ptb(i) or ptb(i-1) ) [un>0 or <0] 
    70       !!               zupsv = e1v*e3v  vn  (tb(j) or tb(j-1) ) [vn>0 or <0] 
    71       !!       * mixed upstream / centered horizontal advection scheme 
     62      !!         Part I : advective fluxes 
     63      !!       * centered flux:   zcenu = pun  mi(ptn)         (idem for v and w) 
     64      !!       * upstream flux:   zupsu = pun  (ptb(i) or ptb(i-1) ) [un>0 or <0] 
     65      !!       * mixed upstream / centered flux: 
    7266      !!               zcofi = max(zind(i+1), zind(i)) 
    73       !!               zcofj = max(zind(j+1), zind(j)) 
    7467      !!               zwx = zcofi * zupsu + (1-zcofi) * zcenu 
    75       !!               zwy = zcofj * zupsv + (1-zcofj) * zcenv 
     68      !!         Part III : advective trend 
    7669      !!       * horizontal advective trend (divergence of the fluxes) 
    77       !!               zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] } 
    78       !!       * Add this trend now to the general trend of tracer (pta): 
    79       !!               pta = pta + zta  
    80       !!       * trend diagnostic (lk_trdtra=T): the trend is 
    81       !!      saved for diagnostics. The trends saved is expressed as 
    82       !!      Uh.gradh(T), i.e. 
    83       !!                     save trend = zta + ptn divn 
    84       !!         In addition, the advective trend in the two horizontal direc- 
    85       !!      tion is also re-computed as Uh gradh(T). Indeed hadt+ptn divn is 
    86       !!      equal to (in s-coordinates, and similarly in z-coord.): 
    87       !!         zta+ptn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[ptn] ) 
    88       !!                                       +mj-1( e1v*e3v  vn  mj[ptn] )  } 
    89       !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so 
    90       !!      they vanish from the expression of the flux and divergence. 
    91       !! 
    92       !!         Part II : vertical advection 
    93       !!      the advective trend is computed as follows : 
    94       !!            zta = 1/e3t dk+1[ zwx ] 
    95       !!      where the vertical advective flux, zwx, is given by : 
    96       !!            zwx = zcofk * zupst + (1-zcofk) * zcent 
    97       !!      with 
    98       !!        zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0] 
    99       !!        zcenu = centered flux = wn * mk(ptn) 
    100       !!         The surface boundary condition is : 
    101       !!      rigid-lid (lk_dynspg_frd = T) : zero advective flux 
    102       !!      free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * ptn(:,:,1) 
    103       !!         Add this trend now to the general trend of tracer (ta,sa): 
    104       !!            (ta,sa) = (ta,sa) + ( zta , zsa ) 
    105       !!         Trend diagnostic ('key_trdtra' defined): the trend is 
    106       !!      saved for diagnostics. The trends saved is expressed as : 
    107       !!             save trend =  w.gradz(T) = zta - ptn divn. 
    108       !! 
    109       !! ** Action :  - update (ta,sa) with the now advective tracer trends 
     70      !!         added to the general trend of tracer (pta): 
     71      !!               pta = pta - 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] + dk[zwz] } 
     72      !!       * trend diagnostic (lk_trdtra=T): the advective fluxes are 
     73      !!      where the 3 part of advective trend are recomputed as U.grad[T] 
     74      !!      and then used in on-line diagnostics. 
     75      !! 
     76      !!         NB: The surface boundary condition is no flux for both 
     77      !!      rigid-lid (lk_dynspg_frd = T) and non-linear free-surface 
     78      !!      but a flux (pwn(1) * ptn(1)) cross the surface in linear 
     79      !!      free-surface case (lk_dynspg_fsc = T) 
     80      !! 
     81      !! ** Action :  - update pta with the 2nd order centered advective trends 
    11082      !!              - trend diagnostics (lk_trdtra=T) 
    11183      !!---------------------------------------------------------------------- 
     
    11890      !! 
    11991      INTEGER  ::   ji, jj, jk                                    ! dummy loop indices 
    120       REAL(wp) ::   zbtr, zta, zhw, ze3tr                         ! temporary scalars 
    121       REAL(wp) ::   zcofi, zfui, zcenut, zupsut, zfp_ui, zfm_ui   !    "         " 
    122       REAL(wp) ::   zcofj, zfvj, zcenvt, zupsvt, zfp_vj, zfm_vj   !    "         " 
    123       REAL(wp) ::   zcofk,       zcent , zupst , zfp_w , zfm_w    !    "         " 
    124       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zind   ! 3D workspace  
     92      REAL(wp) ::   zbtr                                          ! temporary scalars 
     93      REAL(wp) ::   zcofi, zcenut, zupsut, zfp_ui, zfm_ui         !    "         " 
     94      REAL(wp) ::   zcofj, zcenvt, zupsvt, zfp_vj, zfm_vj         !    "         " 
     95      REAL(wp) ::   zcofk, zcent , zupst , zfp_w , zfm_w          !    "         " 
     96      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy, zwz, zind   ! 3D workspace  
    12597      !!---------------------------------------------------------------------- 
    12698 
     
    129101         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme' 
    130102         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    131          !  
    132          btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
    133103      ENDIF 
    134104 
     
    150120      END DO 
    151121 
    152       ! I. Horizontal advection 
    153       ! ----------------------- 
    154  
    155       !                                                ! =============== 
    156       DO jk = 1, jpkm1                                 ! Horizontal slab 
    157          !                                             ! =============== 
    158          ! 
    159          DO jj = 1, jpjm1                              !  Horizontal advective fluxes 
     122      ! I. Advective fluxes 
     123      ! ------------------- 
     124      DO jk = 1, jpkm1                                 ! Horizontal advective fluxes 
     125         DO jj = 1, jpjm1 
    160126            DO ji = 1, fs_jpim1   ! vector opt. 
    161127               ! upstream indicator 
    162128               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) ) 
    163129               zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) ) 
    164                ! volume fluxes * 1/2 
    165 #if defined key_zco 
    166                zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk) 
    167                zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk) 
    168 #else 
    169                zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 
    170                zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 
    171 #endif 
    172130               ! upstream scheme 
    173                zfp_ui = zfui + ABS( zfui ) 
    174                zfp_vj = zfvj + ABS( zfvj ) 
    175                zfm_ui = zfui - ABS( zfui ) 
    176                zfm_vj = zfvj - ABS( zfvj ) 
     131               zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 
     132               zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 
     133               zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 
     134               zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 
    177135               zupsut = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk) 
    178136               zupsvt = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk) 
    179137               ! centered scheme 
    180                zcenut = zfui * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) ) 
    181                zcenvt = zfvj * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) ) 
     138               zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) ) 
     139               zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) ) 
    182140               ! mixed centered / upstream scheme 
    183                zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut 
    184                zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt 
    185             END DO 
    186          END DO 
    187  
    188          DO jj = 2, jpjm1                              !  horizontal tracer flux divergence added to the general trend 
    189             DO ji = fs_2, fs_jpim1   ! vector opt. 
    190 #if defined key_zco 
    191                zbtr = btr2(ji,jj) 
    192 #else 
    193                zbtr = btr2(ji,jj) / fse3t(ji,jj,jk) 
    194 #endif 
    195                ! horizontal advective trends 
    196                zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
    197                   &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
    198                ! add it to the general tracer trends 
    199                pta(ji,jj,jk) = pta(ji,jj,jk) + zta 
    200             END DO 
    201          END DO 
    202          !                                             ! =============== 
    203       END DO                                           !   End of slab 
    204       !                                                ! =============== 
    205  
    206       !  Save the horizontal advective trends for diagnostic 
    207       ! ----------------------------------------------------- 
    208       IF( l_trdtra ) THEN 
    209          CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn ) 
    210          CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn ) 
    211       ENDIF 
    212  
    213       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - had: ', mask1=tmask, clinfo3=cdtype ) 
    214  
    215  
    216       ! "Poleward" heat and salt transport  
    217       ! ---------------------------------- 
    218       IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
    219          IF( lk_zco ) THEN 
    220             DO jk = 1, jpkm1 
    221                DO jj = 2, jpjm1 
    222                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    223                     zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 
    224                   END DO 
    225                END DO 
    226             END DO 
    227          ENDIF 
    228          IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
    229          IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) ) 
    230       ENDIF 
    231  
    232  
    233       ! II. Vertical advection 
    234       ! ---------------------- 
    235  
    236       IF( lk_dynspg_rl .OR. lk_vvl ) THEN         ! rigid lid or non-linear free surface 
    237          zwx(:,:, 1 ) = 0.e0                           ! Surface value : zero flux 
    238          zwx(:,:,jpk) = 0.e0                           ! Bottom  value : flux set to zero 
    239       ELSE                                        ! linear free surface 
    240          zwx(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)        ! Surface :  : advection through z=0 
    241          zwx(:,:,jpk) = 0.e0                           ! Bottom  : flux set to zero 
    242       ENDIF 
    243  
    244       DO jk = 2, jpk                                   ! Vertical advective fluxes (at w-point) 
     141               zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut ) 
     142               zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt ) 
     143            END DO 
     144         END DO 
     145      END DO 
     146      !                                                ! Vertical advective fluxes 
     147      zwz(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero 
     148      !                                                     ! Surface value :  
     149      IF( lk_dynspg_rl .OR. lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                      ! rigid lid or non-linear fs 
     150      ELSE                                  ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)   ! linear free surface  
     151      ENDIF 
     152      ! 
     153      DO jk = 2, jpkm1                                      ! interior values 
    245154         DO jj = 2, jpjm1 
    246155            DO ji = fs_2, fs_jpim1   ! vector opt. 
    247156               ! upstream indicator 
    248157               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
    249                ! velocity * 1/2 
    250                zhw = 0.5 * pwn(ji,jj,jk) 
    251158               ! upstream scheme 
    252                zfp_w = zhw + ABS( zhw ) 
    253                zfm_w = zhw - ABS( zhw ) 
     159               zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     160               zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    254161               zupst = zfp_w * ptb(ji,jj,jk) + zfm_w * ptb(ji,jj,jk-1) 
    255162               ! centered scheme 
    256                zcent = zhw * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) 
     163               zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) 
    257164               ! mixed centered / upstream scheme 
    258                zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent 
    259             END DO 
    260          END DO 
    261       END DO 
    262  
    263       DO jk = 1, jpkm1                                 ! Tracer flux divergence at t-point added to the general trend 
    264          DO jj = 2, jpjm1 
     165               zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent ) 
     166            END DO 
     167         END DO 
     168      END DO 
     169 
     170      ! II. Divergence of advective fluxes 
     171      ! ---------------------------------- 
     172      DO jk = 1, jpkm1                    
     173         DO jj = 2, jpjm1       
    265174            DO ji = fs_2, fs_jpim1   ! vector opt. 
    266                ze3tr = 1. / fse3t(ji,jj,jk) 
    267                ! vertical advective trends 
    268                zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
    269                ! add it to the general tracer trends 
    270                pta(ji,jj,jk) =  pta(ji,jj,jk) + zta 
    271             END DO 
    272          END DO 
    273       END DO 
    274  
    275       ! 3. Save the vertical advective trends for diagnostic 
    276       ! ---------------------------------------------------- 
    277       IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptn ) 
    278  
    279       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - zad : ', mask1=tmask, clinfo3=cdtype ) 
     175               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     176               ! advective trends added to the general tracer trends 
     177               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     178                  &                                    + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     179                  &                                    + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  ) 
     180            END DO 
     181         END DO 
     182      END DO 
     183 
     184 
     185      ! III. Diagnostics and control print 
     186      ! ---------------------------------- 
     187      IF( l_trdtra ) THEN                            ! trend diagnostics 
     188         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn ) 
     189         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn ) 
     190         CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwz, pwn, ptn ) 
     191      ENDIF 
     192      !                                              ! "Poleward" heat and salt transport  
     193      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
     194         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
     195         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) ) 
     196      ENDIF 
     197      !                                              ! control print 
     198      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' cen2 - adv: ', mask1=tmask, clinfo3=cdtype ) 
    280199      ! 
    281200   END SUBROUTINE tra_adv_cen2 
  • branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_eiv.F90

    r719 r791  
    55   !!====================================================================== 
    66   !! History :  9.0  !  05-11  (G. Madec)  Original code, from traldf and zdf _iso 
     7   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport 
    78   !!---------------------------------------------------------------------- 
    89#if defined key_traldf_eiv   ||   defined key_esopa 
     
    3233#  include "vectopt_loop_substitute.h90" 
    3334   !!---------------------------------------------------------------------- 
    34    !!  OPA 9.0 , LOCEAN-IPSL (2006)  
    35    !! $Header$  
     35   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)  
     36   !! $Id:$  
    3637   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3738   !!---------------------------------------------------------------------- 
     
    4344      !!                  ***  ROUTINE tra_adv_eiv  *** 
    4445      !!  
    45       !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive  
    46       !!      trend and add it to the general trend of tracer equation. 
    47       !! 
    48       !! ** Method  :   The eddy induced advection is computed from the slope 
     46      !! ** Purpose :   Compute the eddy induced transport and add it to the 
     47      !!              effective transport 
     48      !! 
     49      !! ** Method  :   The eddy induced transport is computed from the slope 
    4950      !!      of iso-neutral surfaces computed in routine ldf_slp as follows: 
    50       !!         zu_eiv =  1/(e2u e3u)   dk[ aeiu e2u mi(wslpi) ] 
    51       !!         zv_eiv =  1/(e1v e3v)   dk[ aeiv e1v mj(wslpj) 
    52       !!         zw_eiv = -1/(e1t e2t) { di[ aeiu e2u mi(wslpi) ] 
    53       !!                               + dj[ aeiv e1v mj(wslpj) ] } 
     51      !!         zu_eiv =   dk[ aeiu e2u mi(wslpi) ] 
     52      !!         zv_eiv =   dk[ aeiv e1v mj(wslpj) ] 
     53      !!         zw_eiv = - { di[ aeiu e2u mi(wslpi) ] + dj[ aeiv e1v mj(wslpj) ] } 
    5454      !!      add the eiv component to the model velocity: 
    5555      !!         p.n = p.n + z._eiv 
     56      !! CAUTION : the horizontal transports not updated along jpi column and jpj row 
     57      !!           the vertical   transports not updated along 1 & jpi columns and 1 & jpj rows 
    5658      !! 
    5759      !! ** Action  : - add to p.n the eiv component 
    5860      !!---------------------------------------------------------------------- 
    5961      INTEGER , INTENT(in   )                         ::   kt    ! ocean time-step index 
    60       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun   ! in : 3 ocean velocity components  
    61       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn   ! out: 3 ocean velocity components 
     62      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun   ! in : 3 ocean transport components  
     63      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn   ! out: 3 ocean transport components 
    6264      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pwn   !      increased by the eiv 
    6365      !! 
     
    8890               zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1) 
    8991 
    90                zu_eiv = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) / fse3u(ji,jj,jk) 
    91                zv_eiv = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) / fse3v(ji,jj,jk) 
     92               zu_eiv = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) 
     93               zv_eiv = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) 
    9294    
    93                pun(ji,jj,jk) = pun(ji,jj,jk) + zu_eiv 
    94                pvn(ji,jj,jk) = pvn(ji,jj,jk) + zv_eiv 
    95 # if defined key_diaeiv 
    96                u_eiv(ji,jj,jk) = zu_eiv 
    97                v_eiv(ji,jj,jk) = zv_eiv 
     95               pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv 
     96               pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv 
     97# if defined key_diaeiv 
     98               u_eiv(ji,jj,jk) = zu_eiv / fse3u(ji,jj,jk) 
     99               v_eiv(ji,jj,jk) = zv_eiv / fse3v(ji,jj,jk) 
    98100# endif 
    99101            END DO 
     
    115117                  zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk) 
    116118 
    117                   zw_eiv = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )  & 
    118                      &                                                / ( e1t(ji,jj)*e2t(ji,jj) ) 
     119                  zw_eiv = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) 
    119120# endif 
    120121                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv 
    121122 
    122123# if defined key_diaeiv 
    123                   w_eiv(ji,jj,jk) = zw_eiv 
     124                  w_eiv(ji,jj,jk) = zw_eiv / ( e1t(ji,jj)*e2t(ji,jj) ) 
    124125# endif 
    125126               END DO 
     
    130131      !                                             ! ================= 
    131132   END SUBROUTINE tra_adv_eiv 
     133 
     134 
     135!!gm   test tra_adv_eiv better (faster) coded?  to be verified 
     136 
     137   SUBROUTINE tra_adv_eiv2( kt, pun, pvn, pwn ) 
     138      !!---------------------------------------------------------------------- 
     139      !!                  ***  ROUTINE tra_adv_eiv  *** 
     140      !!  
     141      !! ** Purpose :   Compute the eddy induced transport and add it to the 
     142      !!              effective transport 
     143      !! 
     144      !! ** Method  :   The eddy induced transport is computed from the slope 
     145      !!              of iso-neutral surfaces (see ldfslp.F90) as follows: 
     146      !!                   zu_eiv =   dk[ aeiu e2u mi(wslpi) ] 
     147      !!                   zv_eiv =   dk[ aeiv e1v mj(wslpj) ] 
     148      !!                   zw_eiv = - { di[ aeiu e2u mi(wslpi) ] + dj[ aeiv e1v mj(wslpj) ] } 
     149      !!              add the eiv component to the model velocity: 
     150      !!                   p.n = p.n + z._eiv 
     151      !! CAUTION : the horizontal transports not updated along jpi column and jpj row 
     152      !!           the vertical   transports not updated along 1 & jpi columns and 1 & jpj rows 
     153      !! 
     154      !! ** Action  : - add to p.n the eiv transport component 
     155      !!---------------------------------------------------------------------- 
     156      INTEGER , INTENT(in   )                         ::   kt    ! ocean time-step index 
     157      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pun   ! in : 3 ocean transport components  
     158      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pvn   ! out: 3 ocean transport components 
     159      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pwn   !      increased by the eiv 
     160      !! 
     161      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
     162      REAL(wp) ::   zuwslpi, zvwslpj           ! temporary scalar 
     163      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsfu, zsfv   ! eiv stream-function in u and v directions 
     164      !!---------------------------------------------------------------------- 
     165 
     166      IF( kt == nit000 ) THEN 
     167         IF(lwp) WRITE(numout,*) 
     168         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection :' 
     169         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component' 
     170      ENDIF 
     171 
     172      ! eiv stream-function in u- and v-directions 
     173      ! NB: UW-point mask at level k is umask(:,:,k) idem form VW-point mask 
     174      zsfu(:,:, 1 ) = 0.e0   ;   zsfv(:,:, 1 ) = 0.e0      ! surface value set to zero 
     175      zsfu(:,:,jpk) = 0.e0   ;   zsfv(:,:,jpk) = 0.e0      ! bottom  value set to zero 
     176      DO jk = 2, jpkm1 
     177         DO jj = 1, jpjm1 
     178            DO ji = 1, fs_jpim1   ! vector opt. 
     179               zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 
     180               zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )  
     181               zsfu(ji,jj,jk) = zuwslpi * e2u(ji,jj) * fsaeiu(ji,jj,jk) * umask(ji,jj,jk) 
     182               zsfv(ji,jj,jk) = zvwslpj * e1v(ji,jj) * fsaeiv(ji,jj,jk) * vmask(ji,jj,jk) 
     183            END DO 
     184         END DO 
     185      END DO 
     186# if defined key_diaeiv 
     187      ! save eiv stream function in the output 
     188!!gm  to be done, u_sfeiv and v_sfeiv not defined    ==> new IOM.... 
     189!!gm  and zsfu, zfv notdefined for jpi column and jpj row 
     190!     u_sfeiv(:,:,:) = zsfu(:,:,:) 
     191!     v_sfeiv(:,:,:) = zsfu(:,:,:) 
     192# endif 
     193 
     194      ! increase the transport with the eiv transport 
     195      DO jk = 1, jpkm1 
     196         DO jj = 1, jpjm1 
     197            DO ji = 1, fs_jpim1   ! vector opt. 
     198               pun(ji,jj,jk) = pun(ji,jj,jk) + ( zsfu(ji,jj,jk) - zsfu(ji,jj,jk+1) ) 
     199               pvn(ji,jj,jk) = pvn(ji,jj,jk) + ( zsfv(ji,jj,jk) - zsfv(ji,jj,jk+1) ) 
     200            END DO 
     201         END DO 
     202      END DO 
     203      DO jk = 2, jpkm1 
     204         DO jj = 2, jpjm1 
     205            DO ji = fs_2, fs_jpim1   ! vector opt. 
     206               pwn(ji,jj,jk) = pwn(ji,jj,jk) - (  zsfu(ji,jj,jk) - zsfu(ji-1,jj  ,jk)                        & 
     207                  &                             + zsfv(ji,jj,jk) - zsfv(ji  ,jj-1,jk)  ) * tmask(ji,jj,jk) 
     208            END DO 
     209         END DO 
     210      END DO 
     211      ! 
     212   END SUBROUTINE tra_adv_eiv2 
    132213 
    133214#else 
  • branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl.F90

    r786 r791  
    44   !! Ocean active tracers:  horizontal & vertical advective trend 
    55   !!====================================================================== 
    6    !! History :       !  06-00  (A.Estublier)  for passive tracers 
    7    !!                 !  01-08  (E.Durand, G.Madec)  adapted for T & S 
    8    !!   NEMO     1.0  !  02-06  (G. Madec)  F90: Free form and module 
    9    !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC 
     6   !! History :  OPA  !  2006-00  (A.Estublier)  for passive tracers 
     7   !!            8.2  !  2001-08  (E.Durand, G.Madec)  adapted for T & S 
     8   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module 
     9   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport 
    1010   !!---------------------------------------------------------------------- 
    1111 
     
    2828   PRIVATE 
    2929 
    30    PUBLIC   tra_adv_muscl  ! routine called by step.F90 
     30   PUBLIC   tra_adv_muscl   ! routine called by step.F90 
    3131 
    3232   !! * Substitutions 
     
    3535   !!---------------------------------------------------------------------- 
    3636   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)  
    37    !! $Id:$  
     37   !! $Id$  
    3838   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3939   !!---------------------------------------------------------------------- 
     
    6666      !! 
    6767      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    68       REAL(wp) ::   zu, zv, zw, zeu, zev   
    69       REAL(wp) ::   zew, zbtr, z2, zstep  
    70       REAL(wp) ::   z0u, z0v, z0w  
    71       REAL(wp) ::   zzwx, zzwy, zalpha  
    72       REAL(wp) ::   zta 
     68      REAL(wp) ::   zu, z0u, zzwx 
     69      REAL(wp) ::   zv, z0v, zzwy  
     70      REAL(wp) ::   zw, z0w  
     71      REAL(wp) ::   zbtr, z2, zdt, zalpha  
    7372      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zwx, zwy, zslpx, zslpy   ! 3D workspace 
    7473      !!---------------------------------------------------------------------- 
     
    8079      ENDIF 
    8180 
    82       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1. 
    83       ELSE                                        ;    z2 = 2. 
     81      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.      ! euler     time-stepping 
     82      ELSE                                        ;    z2 = 2.      ! leap-frog time-stepping 
    8483      ENDIF 
    8584 
    8685      ! I. Horizontal advective fluxes 
    8786      ! ------------------------------ 
    88       ! first guess of the slopes 
    89       ! interior values 
    90       DO jk = 1, jpkm1 
     87      !                                             !-- first guess of the slopes 
     88      zwx(:,:,jpk) = 0.e0   ;   zwy(:,:,jpk) = 0.e0        ! bottom values 
     89      DO jk = 1, jpkm1                                     ! interior values 
    9190         DO jj = 1, jpjm1       
    9291            DO ji = 1, fs_jpim1   ! vector opt. 
     
    9695         END DO 
    9796      END DO 
    98       ! bottom values 
    99       zwx(:,:,jpk) = 0.e0    ;    zwy(:,:,jpk) = 0.e0 
    100  
    101       ! lateral boundary conditions on zwx, zwy   (changed sign) 
    102       CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. ) 
    103  
    104       ! Slopes 
    105       ! interior values 
    106       DO jk = 1, jpkm1 
    107          DO jj = 2, jpj 
    108             DO ji = fs_2, jpi   ! vector opt. 
     97!!gm optimisation (lbc to be suppressed) 
     98      CALL lbc_lnk( zwx, 'U', -1. )                        ! lateral boundary conditions on zwx, zwy   (changed sign) 
     99      CALL lbc_lnk( zwy, 'V', -1. ) 
     100!!gm 
     101 
     102      !                                             !-- Slopes of tracer 
     103      zslpx(:,:,jpk) = 0.e0   ;   zslpy(:,:,jpk) = 0.e0    ! bottom values 
     104      DO jk = 1, jpkm1                                     ! interior values 
     105         DO jj = 2, jpjm1 
     106            DO ji = fs_2, jpim1   ! vector opt. 
    109107               zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   & 
    110                   &           * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) ) 
     108                  &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) ) 
    111109               zslpy(ji,jj,jk) =                    ( zwy(ji,jj,jk) + zwy(ji  ,jj-1,jk) )   & 
    112                   &           * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) ) 
    113             END DO 
    114          END DO 
    115       END DO 
    116       ! bottom values 
    117       zslpx(:,:,jpk) = 0.e0    ;    zslpy(:,:,jpk) = 0.e0 
    118  
    119       ! Slopes limitation 
    120       DO jk = 1, jpkm1 
    121          DO jj = 2, jpj 
    122             DO ji = fs_2, jpi   ! vector opt. 
    123                zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) )   & 
    124                   &            * MIN(    ABS( zslpx(ji  ,jj,jk) ),   & 
    125                   &                   2.*ABS( zwx (ji-1,jj,jk) ),   & 
    126                   &                   2.*ABS( zwx (ji  ,jj,jk) ) ) 
    127                zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) )   & 
    128                   &            * MIN(    ABS( zslpy(ji,jj  ,jk) ),   & 
    129                   &                   2.*ABS( zwy (ji,jj-1,jk) ),   & 
    130                   &                   2.*ABS( zwy (ji,jj  ,jk) ) ) 
     110                  &            * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji  ,jj-1,jk) ) ) 
     111            END DO 
     112         END DO 
     113      END DO 
     114!!gm  :merge the 2 loop (above & below) 
     115      DO jk = 1, jpkm1                                     ! Slopes limitation 
     116         DO jj = 2, jpjm1 
     117            DO ji = fs_2, jpim1   ! vector opt. 
     118               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   & 
     119                  &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   & 
     120                  &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) ) 
     121               zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   & 
     122                  &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   & 
     123                  &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) ) 
    131124            END DO 
    132125         END DO 
    133126      END DO         
    134  
    135       ! Advection terms 
    136       ! interior values 
    137       DO jk = 1, jpkm1 
    138          zstep  = z2 * rdttra(jk) 
    139          DO jj = 2, jpjm1       
    140             DO ji = fs_2, fs_jpim1   ! vector opt. 
    141                ! volume fluxes 
    142 #if defined key_zco 
    143                zeu = e2u(ji,jj)                   * pun(ji,jj,jk) 
    144                zev = e1v(ji,jj)                   * pvn(ji,jj,jk) 
    145 #else 
    146                zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 
    147                zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 
    148 #endif 
    149                ! MUSCL fluxes 
     127!!gm optimisation:  call lbclnk just 1 time here on zslpx & zslpy 
     128!!    CALL lbc_lnk( zslpx, 'U', -1. )                      ! lateral boundary conditions on zwx, zwy   (changed sign) 
     129!!    CALL lbc_lnk( zslpy, 'V', -1. ) 
     130!! 
     131!! and loop below from 1 to jpjm1 and to jpim1 
     132!!gm 
     133 
     134      !                                             !-- MUSCL horizontal advective fluxes 
     135      DO jk = 1, jpkm1                                     ! interior values 
     136         zdt  = z2 * rdttra(jk) 
     137         DO jj = 2, jpjm1       
     138            DO ji = fs_2, fs_jpim1   ! vector opt. 
    150139               z0u = SIGN( 0.5, pun(ji,jj,jk) )             
    151140               zalpha = 0.5 - z0u 
    152                zu  = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj) 
     141               zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 
    153142               zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk) 
    154143               zzwy = ptb(ji  ,jj,jk) + zu * zslpx(ji  ,jj,jk) 
    155                zwx(ji,jj,jk) = zeu * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     144               zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    156145               ! 
    157146               z0v = SIGN( 0.5, pvn(ji,jj,jk) )             
    158147               zalpha = 0.5 - z0v 
    159                zv  = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj) 
     148               zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 
    160149               zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk) 
    161150               zzwy = ptb(ji,jj  ,jk) + zv * zslpy(ji,jj  ,jk) 
    162                zwy(ji,jj,jk) = zev * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    163             END DO 
    164          END DO 
    165       END DO 
    166  
    167 !!gm bug?   there is too many lbc: this have to be checked 
    168       ! lateral boundary conditions on zwx, zwy   (changed sign) 
     151               zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     152            END DO 
     153         END DO 
     154      END DO 
     155!!gm optimisation (lbc to be suppressed) 
     156      !                                                    ! lateral boundary conditions on zwx, zwy   (changed sign) 
    169157      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. ) 
    170  
     158!!gm 
     159 
     160      !                                             !-- MUSCL advective trend 
    171161      ! Tracer flux divergence at t-point added to the general trend 
    172162      DO jk = 1, jpkm1 
    173163         DO jj = 2, jpjm1       
    174164            DO ji = fs_2, fs_jpim1   ! vector opt. 
    175 #if defined key_zco 
    176                zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 
    177 #else 
    178165               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    179 #endif 
    180                ! horizontal advective trends 
    181                zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    182                   &           + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )  
    183                ! add it to the general tracer trends 
    184                pta(ji,jj,jk) = pta(ji,jj,jk) + zta 
     166               ! horizontal advective trends added to the general tracer trends 
     167               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     168                  &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )  
    185169            END DO 
    186170        END DO 
    187171      END DO         
    188172 
    189       IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype ) 
    190  
    191       ! Save the horizontal advective trends for diagnostics 
    192       IF( l_trdtra ) THEN 
     173      IF( l_trdtra ) THEN                           !-- trend diagnostics 
    193174         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptb ) 
    194175         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptb ) 
    195176      ENDIF 
    196  
     177      !                                             !-- "Poleward" heat or salt transports 
    197178      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
    198          IF( lk_zco ) THEN 
    199             DO jk = 1, jpkm1 
    200                DO jj = 2, jpjm1 
    201                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    202                     zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 
    203                   END DO 
    204                END DO 
    205             END DO 
    206          ENDIF 
    207179         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
    208180         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) ) 
    209181      ENDIF 
     182      !                                             !-- control print 
     183      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype ) 
    210184 
    211185 
     
    213187      ! ----------------------------- 
    214188       
    215       ! First guess of the slope 
    216       ! interior values 
    217       DO jk = 2, jpkm1 
     189      !                                             !-- first guess of the slopes 
     190      zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0    ! surface & bottom boundary conditions 
     191      DO jk = 2, jpkm1                                     ! interior values 
    218192         zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1) - ptb(:,:,jk) ) 
    219193      END DO 
    220       ! surface & bottom boundary conditions 
    221       zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0 
    222  
    223       ! Slopes 
    224       DO jk = 2, jpkm1 
     194 
     195      !                                             !-- Slopes of tracer 
     196      zslpx(:,:,1) = 0.e0                                  ! surface values 
     197      DO jk = 2, jpkm1                                     ! interior value 
    225198         DO jj = 1, jpj 
    226199            DO ji = 1, jpi 
     
    230203         END DO 
    231204      END DO 
    232  
    233       ! Slopes limitation 
    234       DO jk = 2, jpkm1        ! interior values 
     205!!gm : merge the above and below loops 
     206      !                                             !-- Slopes limitation 
     207      DO jk = 2, jpkm1                                     ! interior values 
    235208         DO jj = 1, jpj 
    236209            DO ji = 1, jpi 
    237                zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) )   & 
    238                   &            * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
    239                   &                   2.*ABS( zwx (ji,jj,jk+1) ),   & 
    240                   &                   2.*ABS( zwx (ji,jj,jk  ) ) ) 
    241             END DO 
    242          END DO 
    243       END DO 
    244       zslpx(:,:,1) = 0.e0      ! surface values 
    245  
    246       ! vertical advective flux 
    247       DO jk = 1, jpkm1        ! interior values 
    248          zstep  = z2 * rdttra(jk) 
    249          DO jj = 2, jpjm1       
    250             DO ji = fs_2, fs_jpim1   ! vector opt. 
    251                zew = pwn(ji,jj,jk+1) 
     210               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
     211                  &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   & 
     212                  &                                                 2.*ABS( zwx  (ji,jj,jk  ) )  ) 
     213            END DO 
     214         END DO 
     215      END DO 
     216 
     217      !                                             !-- vertical advective flux 
     218      !                                                    ! surface values  (bottom already set to zero) 
     219      IF( lk_dynspg_rl .OR. lk_vvl) THEN    ;   zwx(:,:, 1 ) = 0.e0                      ! rigid lid or variable volume 
     220      ELSE                                  ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)   ! linear free surface 
     221      ENDIF 
     222      DO jk = 1, jpkm1                                     ! interior values 
     223         zdt  = z2 * rdttra(jk) 
     224         DO jj = 2, jpjm1       
     225            DO ji = fs_2, fs_jpim1   ! vector opt. 
    252226               z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
    253227               zalpha = 0.5 + z0w 
    254                zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1) 
     228               zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt / (e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 
    255229               zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1) 
    256230               zzwy = ptb(ji,jj,jk  ) + zw * zslpx(ji,jj,jk  ) 
    257                zwx(ji,jj,jk+1) = zew * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    258             END DO 
    259          END DO 
    260       END DO 
    261       !                       ! surface values 
    262       IF( lk_dynspg_rl .OR. lk_vvl) THEN                ! rigid lid or variable volume: flux set to zero 
    263          zwx(:,:, 1 ) = 0.e0                                 ! surface  
    264          zwx(:,:,jpk) = 0.e0                                 ! bottom  
    265       ELSE                                              ! free surface 
    266          zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)              ! Surface 
    267          zwx(:,:,jpk) = 0.e0                                 ! bottom  
    268  
    269       ENDIF 
    270  
    271       ! Compute & add the vertical advective trend 
     231               zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     232            END DO 
     233         END DO 
     234      END DO 
     235 
     236      !                                             !-- vertical advective trend 
    272237      DO jk = 1, jpkm1 
    273238         DO jj = 2, jpjm1       
    274239            DO ji = fs_2, fs_jpim1   ! vector opt. 
    275                zbtr = 1. / fse3t(ji,jj,jk) 
    276                ! horizontal advective trends 
    277                zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
    278                ! add it to the general tracer trends 
    279                pta(ji,jj,jk) =  pta(ji,jj,jk) + zta 
    280             END DO 
    281          END DO 
    282       END DO 
    283  
    284       ! Save the vertical advective trends for diagnostic 
    285       ! ------------------------------------------------- 
     240               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     241               ! horizontal advective trends added to the general tracer trends 
     242               pta(ji,jj,jk) =  pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
     243            END DO 
     244         END DO 
     245      END DO 
     246 
     247      !                                             !-- trend diagnostic 
    286248      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptb ) 
    287  
     249      !                                             !-- control print 
    288250      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - zad: ', mask1=tmask, clinfo3=cdtype ) 
    289251      ! 
  • branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90

    r786 r791  
    44   !! Ocean active tracers:  horizontal & vertical advective trend 
    55   !!============================================================================== 
    6    !! History :  1.0  !  02-06  (G. Madec) from traadv_muscl 
    7    !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC 
     6   !! History :  1.0  !  2002-06  (G. Madec) from traadv_muscl 
     7   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport 
    88   !!---------------------------------------------------------------------- 
    99 
     
    2626   PRIVATE 
    2727 
    28    !! * Accessibility 
    2928   PUBLIC tra_adv_muscl2        ! routine called by step.F90 
    3029 
     
    3433   !!---------------------------------------------------------------------- 
    3534   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)  
    36    !! $Id:$  
     35   !! $Id$  
    3736   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3837   !!---------------------------------------------------------------------- 
     
    6665      !! 
    6766      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    68       REAL(wp) ::   zu, zv, zw, zeu, zev 
    69       REAL(wp) ::   zew, zbtr, zstep, z2  
     67      REAL(wp) ::   zu, zv, zw 
     68      REAL(wp) ::   zbtr, zstep, z2  
    7069      REAL(wp) ::   z0u, z0v, z0w 
    71       REAL(wp) ::   zzwx, zzwy, zalpha, zta 
     70      REAL(wp) ::   zzwx, zzwy, zalpha 
    7271      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zwx, zwy, zslpx, zslpy   ! 3D workspace 
    7372      !!---------------------------------------------------------------------- 
     
    9897      zwx(:,:,jpk) = 0.e0    ;    zwy(:,:,jpk) = 0.e0        ! bottom values 
    9998 
     99!!gm  this lbc_lnk can be omitted 
    100100      ! lateral boundary conditions on zwx, zwy   (changed sign) 
    101101      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. ) 
     
    117117         DO jj = 2, jpj 
    118118            DO ji = fs_2, jpi   ! vector opt. 
    119                zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) )   & 
    120                   &            * MIN(    ABS( zslpx(ji  ,jj,jk) ),   & 
    121                   &                   2.*ABS( zwx (ji-1,jj,jk) ),   & 
    122                   &                   2.*ABS( zwx (ji  ,jj,jk) ) ) 
    123                zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) )   & 
    124                   &            * MIN(    ABS( zslpy(ji,jj  ,jk) ),   & 
    125                   &                   2.*ABS( zwy (ji,jj-1,jk) ),   & 
    126                   &                   2.*ABS( zwy (ji,jj  ,jk) ) ) 
     119               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   & 
     120                  &                                                 2.*ABS( zwx  (ji-1,jj,jk) ),   & 
     121                  &                                                 2.*ABS( zwx  (ji  ,jj,jk) ) ) 
     122               zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN(    ABS( zslpy(ji,jj  ,jk) ),   & 
     123                  &                                                 2.*ABS( zwy  (ji,jj-1,jk) ),   & 
     124                  &                                                 2.*ABS( zwy  (ji,jj  ,jk) ) ) 
    127125            END DO 
    128126         END DO 
     
    134132         DO jj = 2, jpjm1       
    135133            DO ji = fs_2, fs_jpim1   ! vector opt. 
    136                ! volume fluxes 
    137 #if defined key_zco 
    138                zeu = e2u(ji,jj)                   * pun(ji,jj,jk) 
    139                zev = e1v(ji,jj)                   * pvn(ji,jj,jk) 
    140 #else 
    141                zeu = e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 
    142                zev = e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 
    143 #endif 
    144                ! MUSCL fluxes 
    145                z0u = SIGN( 0.5, pun(ji,jj,jk) )             
     134               z0u  = SIGN( 0.5, pun(ji,jj,jk) )             
    146135               zalpha = 0.5 - z0u 
    147                zu  = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj) 
    148                zzwx = ptb(ji+1,jj,jk) + zu*zslpx(ji+1,jj,jk) 
    149                zzwy = ptb(ji  ,jj,jk) + zu*zslpx(ji  ,jj,jk) 
    150                zwx(ji,jj,jk) = zeu * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     136               zu   = z0u - 0.5 * pun(ji,jj,jk) * zstep / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 
     137               zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk) 
     138               zzwy = ptb(ji  ,jj,jk) + zu * zslpx(ji  ,jj,jk) 
     139               zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    151140               ! 
    152141               z0v = SIGN( 0.5, pvn(ji,jj,jk) )             
    153142               zalpha = 0.5 - z0v 
    154                zv  = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj) 
    155                zzwx = ptb(ji,jj+1,jk) + zv*zslpy(ji,jj+1,jk) 
    156                zzwy = ptb(ji,jj  ,jk) + zv*zslpy(ji,jj  ,jk) 
    157                zwy(ji,jj,jk) = zev * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
     143               zv   = z0v - 0.5 * pvn(ji,jj,jk) * zstep / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 
     144               zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk) 
     145               zzwy = ptb(ji,jj  ,jk) + zv * zslpy(ji,jj  ,jk) 
     146               zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    158147            END DO 
    159148         END DO 
     
    166155        DO jj = 2, jpjm1 
    167156            DO ji = fs_2, fs_jpim1   ! vector opt. 
    168 #if defined key_zco 
    169157               IF( umask(ji,jj,jk) == 0. ) THEN 
    170158                  IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN 
    171                      zwx(ji+1,jj,jk) = e2u(ji+1,jj) * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5 
     159                     zwx(ji+1,jj,jk) = pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5 
    172160                  ENDIF 
    173161                  IF( pun(ji-1,jj,jk) < 0. ) THEN 
    174                      zwx(ji-1,jj,jk) = e2u(ji-1,jj) * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji  ,jj,jk) ) * 0.5 
     162                     zwx(ji-1,jj,jk) = pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji  ,jj,jk) ) * 0.5 
    175163                  ENDIF 
    176164               ENDIF 
    177165               IF( vmask(ji,jj,jk) == 0. ) THEN 
    178166                  IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN 
    179                      zwy(ji,jj+1,jk) = e1v(ji,jj+1) * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5 
     167                     zwy(ji,jj+1,jk) = pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5 
    180168                  ENDIF 
    181169                  IF( pvn(ji,jj-1,jk) < 0. ) THEN 
    182                      zwy(ji,jj-1,jk) = e1v(ji,jj-1) * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji  ,jj,jk) ) * 0.5 
     170                     zwy(ji,jj-1,jk) = pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji  ,jj,jk) ) * 0.5 
    183171                  ENDIF 
    184172               ENDIF 
    185 #else 
    186                IF( umask(ji,jj,jk) == 0. ) THEN 
    187                   IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN 
    188                      zwx(ji+1,jj,jk) = e2u(ji+1,jj)* fse3u(ji+1,jj,jk)   & 
    189                         &            * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk) + ptn(ji+2,jj,jk) ) * 0.5 
    190                   ENDIF 
    191                   IF( pun(ji-1,jj,jk) < 0. ) THEN 
    192                      zwx(ji-1,jj,jk) = e2u(ji-1,jj)* fse3u(ji-1,jj,jk)   & 
    193                         &            * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk) + ptn(ji  ,jj,jk) ) * 0.5 
    194                   ENDIF 
    195                ENDIF 
    196                IF( vmask(ji,jj,jk) == 0. ) THEN 
    197                   IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN 
    198                      zwy(ji,jj+1,jk) = e1v(ji,jj+1) * fse3v(ji,jj+1,jk)   & 
    199                         &            * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk) + ptn(ji,jj+2,jk) ) * 0.5 
    200                   ENDIF 
    201                   IF( pvn(ji,jj-1,jk) < 0. ) THEN 
    202                      zwy(ji,jj-1,jk) = e1v(ji,jj-1)* fse3v(ji,jj-1,jk)   & 
    203                         &            * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk) + ptn(ji  ,jj,jk) ) * 0.5 
    204                   ENDIF 
    205                ENDIF 
    206 #endif 
    207173            END DO 
    208174         END DO 
     
    217183         DO jj = 2, jpjm1       
    218184            DO ji = fs_2, fs_jpim1   ! vector opt. 
    219 #if defined key_zco 
    220                zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) ) 
    221 #else 
    222185               zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 
    223 #endif 
    224                ! horizontal advective trends 
    225                zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    226                   &           + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )  
    227                ! add it to the general tracer trends 
    228                pta(ji,jj,jk) = pta(ji,jj,jk) + zta 
     186               ! horizontal advective trends added to the general tracer trends 
     187               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     188                  &                                   + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )  
    229189            END DO 
    230190        END DO 
     
    270230            DO ji = 1, jpi 
    271231               zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   & 
    272                   &           * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 
     232                  &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 
    273233            END DO 
    274234         END DO 
     
    279239         DO jj = 1, jpj 
    280240            DO ji = 1, jpi 
    281                zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) )   & 
    282                   &           * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
    283                   &                  2.*ABS( zwx (ji,jj,jk+1) ),   & 
    284                   &                  2.*ABS( zwx (ji,jj,jk  ) ) ) 
     241               zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji,jj,jk  ) ),   & 
     242                  &                                                 2.*ABS( zwx  (ji,jj,jk+1) ),   & 
     243                  &                                                 2.*ABS( zwx  (ji,jj,jk  ) ) ) 
    285244            END DO 
    286245         END DO 
     
    293252         DO jj = 2, jpjm1       
    294253            DO ji = fs_2, fs_jpim1   ! vector opt. 
    295                zew = pwn(ji,jj,jk+1) 
    296254               z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
    297255               zalpha = 0.5 + z0w 
    298                zw  = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1) 
    299                zzwx = ptb(ji,jj,jk+1) + zw*zslpx(ji,jj,jk+1) 
    300                zzwy = ptb(ji,jj,jk  ) + zw*zslpx(ji,jj,jk  ) 
    301                zwx(ji,jj,jk+1) = zew * ( zalpha * zzwx + (1.-zalpha)*zzwy ) 
     256               zw  = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / ( e1t(ji,jj) *e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 
     257               zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1) 
     258               zzwy = ptb(ji,jj,jk  ) + zw * zslpx(ji,jj,jk  ) 
     259               zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha)*zzwy ) 
    302260            END DO 
    303261         END DO 
     
    315273      END DO 
    316274 
    317       ! surface values 
    318       IF( lk_dynspg_rl .OR. lk_vvl ) THEN              ! rigid lid or variable volume: flux set to zero 
    319          zwx(:,:, 1 ) = 0.e0                                 ! surface 
    320          zwx(:,:,jpk) = 0.e0                                 ! bottom  
    321       ELSE                                             ! free surface 
    322          zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)               ! surface 
    323          zwx(:,:,jpk) = 0.e0                                 ! bottom  
    324       ENDIF 
     275      zwx(:,:,jpk) = 0.e0                                       ! bottom value 
     276      !                                                         ! surface values 
     277      IF( lk_dynspg_rl .OR. lk_vvl) THEN    ;   zwx(:,:, 1 ) = 0.e0                      ! rigid lid or variable volume 
     278      ELSE                                  ;   zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)   ! linear free surface 
     279      ENDIF  
    325280 
    326281 
     
    330285            DO ji = fs_2, fs_jpim1   ! vector opt. 
    331286               zbtr = 1. / fse3t(ji,jj,jk) 
    332                ! horizontal advective trends 
    333                zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
    334                ! add it to the general tracer trends 
    335                pta(ji,jj,jk) =  pta(ji,jj,jk) + zta 
     287               ! horizontal advective trends added to the general tracer trends 
     288               pta(ji,jj,jk) =  pta(ji,jj,jk) - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
    336289            END DO 
    337290         END DO 
  • branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r786 r791  
    55   !!====================================================================== 
    66   !! History :  1.0  !  06-09  (G. Reffray)  Original code 
    7    !!            2.4  !  08-01  (G. Madec)  Merge TRA-TRC 
     7   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport 
    88   !!---------------------------------------------------------------------- 
    99 
     
    3030   PUBLIC tra_adv_qck    ! routine called by traadv.F90 
    3131 
    32    REAL(wp), DIMENSION(jpi,jpj)     ::   zbtr2 
    3332   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   sl 
    34    REAL(wp) ::   cst1, cst2, dt, coef1                  ! temporary scalars 
    35    INTEGER  ::   ji, jj, jk                             ! dummy loop indices 
    3633 
    3734   !! * Substitutions 
     
    4037   !!---------------------------------------------------------------------- 
    4138   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)  
    42    !! $Id:$  
     39   !! $Id$  
    4340   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4441   !!---------------------------------------------------------------------- 
     
    9491      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend  
    9592      !! 
    96       REAL(wp) :: z2                                   ! temporary scalar  
     93      INTEGER  ::   ji, jj, jk        ! dummy loop indices 
     94      REAL(wp) ::   zdt, z2, zcoef1   ! temporary scalars 
     95      REAL(wp) ::   zbtr              ! temporary scalars 
    9796      !!---------------------------------------------------------------------- 
    9897 
     
    101100         IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3st order quickest advection scheme' 
    102101         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    103  
    104          zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
    105          cst1 = 1./12. 
    106          cst2 = 2./3. 
    107102      ENDIF 
    108103 
    109       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1. 
    110       ELSE                                        ;    z2 = 2. 
     104      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.      ! euler     time-stepping 
     105      ELSE                                        ;    z2 = 2.      ! leap-frog time-stepping 
    111106      ENDIF 
    112107 
     
    116111      !--------------------------------------------------------------- 
    117112 
     113!!gm : optimisation: sl compute twice (for t & s, and even more for trc) 
     114!!gm   note that sl is in permanent memory be used as workspace in the vertical part ! 
    118115      sl(:,:,:) = 100. 
    119116 
    120                                                       ! =============== 
    121       DO jk = 1, jpkm1                                ! Horizontal slab 
    122         !                                             ! =============== 
    123         dt  = z2 * rdttra(jk) 
    124         DO jj = 2, jpjm1 
    125           DO ji = 2, fs_jpim1   ! vector opt. 
    126             coef1 = 1.e-12 
    127             IF (pun(ji-1,jj  ,jk  ).LT.0.) coef1 = coef1 + ABS(pun(ji-1,jj  ,jk  ))*dt/e1t(ji,jj) 
    128             IF (pun(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pun(ji  ,jj  ,jk  ))*dt/e1t(ji,jj) 
    129             IF (pvn(ji  ,jj-1,jk  ).LT.0.) coef1 = coef1 + ABS(pvn(ji  ,jj-1,jk  ))*dt/e2t(ji,jj) 
    130             IF (pvn(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pvn(ji  ,jj  ,jk  ))*dt/e2t(ji,jj) 
    131             IF (pwn(ji  ,jj  ,jk+1).LT.0.) coef1 = coef1 + ABS(pwn(ji  ,jj  ,jk+1))*dt/(fse3t(ji,jj,jk)) 
    132             IF (pwn(ji  ,jj  ,jk  ).GT.0.) coef1 = coef1 + ABS(pwn(ji  ,jj  ,jk  ))*dt/(fse3t(ji,jj,jk)) 
    133             sl(ji,jj,jk) = 1./coef1 
    134             sl(ji,jj,jk) = MIN(100.,sl(ji,jj,jk)) 
    135             sl(ji,jj,jk) = MAX(1.  ,sl(ji,jj,jk)) 
    136           ENDDO 
    137         ENDDO 
    138       ENDDO 
    139       !--- Lateral boundary conditions on the limiter slope 
    140       CALL lbc_lnk(   sl(:,:,:), 'T', 1. ) 
     117      DO jk = 1, jpkm1 
     118         zdt  = z2 * rdttra(jk) 
     119         DO jj = 2, jpjm1 
     120            DO ji = 2, fs_jpim1   ! vector opt. 
     121               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) *fse3t(ji,jj,jk) ) 
     122               zcoef1 = 1.e-12 
     123               IF( pun(ji-1,jj  ,jk  ) < 0.e0 )   zcoef1 = zcoef1 + ABS( pun(ji-1,jj  ,jk  ) ) * zdt * zbtr 
     124               IF( pun(ji  ,jj  ,jk  ) > 0.e0 )   zcoef1 = zcoef1 + ABS( pun(ji  ,jj  ,jk  ) ) * zdt * zbtr 
     125               IF( pvn(ji  ,jj-1,jk  ) < 0.e0 )   zcoef1 = zcoef1 + ABS( pvn(ji  ,jj-1,jk  ) ) * zdt * zbtr 
     126               IF( pvn(ji  ,jj  ,jk  ) > 0.e0 )   zcoef1 = zcoef1 + ABS( pvn(ji  ,jj  ,jk  ) ) * zdt * zbtr 
     127               IF( pwn(ji  ,jj  ,jk+1) < 0.e0 )   zcoef1 = zcoef1 + ABS( pwn(ji  ,jj  ,jk+1) ) * zdt * zbtr 
     128               IF( pwn(ji  ,jj  ,jk  ) > 0.e0 )   zcoef1 = zcoef1 + ABS( pwn(ji  ,jj  ,jk  ) ) * zdt * zbtr 
     129               sl(ji,jj,jk) = 1. / zcoef1 
     130               sl(ji,jj,jk) = MIN( 100., sl(ji,jj,jk) ) 
     131               sl(ji,jj,jk) = MAX( 1.  , sl(ji,jj,jk) ) 
     132            END DO 
     133         END DO 
     134      END DO 
     135      CALL lbc_lnk( sl(:,:,:), 'T', 1. )      ! Lateral boundary conditions on the limiter slope 
    141136 
    142137 
     
    146141      CALL tra_adv_qck_hor( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 ) 
    147142 
     143      !                                       ! control print 
    148144      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - had: ', mask1=tmask, clinfo3=cdtype ) 
    149145 
     
    152148      !------------------------------------------------------------------------- 
    153149 
    154       CALL tra_adv_qck_ver( pwn, ptn , pta, z2 ) 
    155  
     150      CALL tra_adv_qck_ver( pwn, ptn , pta ) 
     151 
     152      !                                       ! control print 
    156153      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' qck - zad: ', mask1=tmask, clinfo3=cdtype ) 
    157154      ! 
     
    159156 
    160157 
    161    SUBROUTINE tra_adv_qck_hor ( kt, cdtype, ktra, pun, pvn, ptb, pta, z2 ) 
    162       !!---------------------------------------------------------------------- 
     158   SUBROUTINE tra_adv_qck_hor ( kt, cdtype, ktra, pun, pvn, ptb, pta, pz2 ) 
     159      !!---------------------------------------------------------------------- 
     160      !!                  ***  ROUTINE tra_adv_qck_hor  *** 
     161      !! 
     162      !! ** Purpose :    
    163163      !! 
    164164      !!---------------------------------------------------------------------- 
     
    166166      CHARACTER(len=3), INTENT(in   )                         ::   cdtype      ! =TRA or TRC (tracer indicator) 
    167167      INTEGER         , INTENT(in   )                         ::   ktra        ! tracer index 
    168       REAL(wp)        , INTENT(in   )                         ::   z2          ! ??? 
     168      REAL(wp)        , INTENT(in   )                         ::   pz2         ! =1 or 2 (euler or leap-frog) 
    169169      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn    ! horizontal effective velocity 
    170170      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb         ! before tracer field 
    171171      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta         ! tracer trend 
    172  
    173       REAL(wp) ::   za, zbtr, e1, e2, c, dir, fu, fc, fd   ! temporary scalars 
    174       REAL(wp) ::   coef2, coef3, fho, mask, dx 
    175       REAL(wp), DIMENSION(jpi,jpj) ::  zee 
    176       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmask, zlap, dwst, lim  
     172      !! 
     173      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     174      REAL(wp) ::   zdt                                     ! temporary scalars 
     175      REAL(wp) ::   zbtr, zc, zdir, zfu, zfc, zfd           ! temporary scalars 
     176      REAL(wp) ::   zcoef1, zcoef2, zcoef3, zfho, zmsk, zdx 
     177      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmask, zlap, zdwst, zlim  
    177178      !!---------------------------------------------------------------------- 
    178179 
     
    182183      zmask(:,:,:)= 0.e0                                                 
    183184      zlap (:,:,:)= 0.e0                                                   
    184       dwst (:,:,:)= 0.e0                                                   
    185       lim (:,:,:)= 0.e0      
     185      zdwst(:,:,:)= 0.e0                                                   
     186      zlim (:,:,:)= 0.e0      
    186187                                             
    187188      !---------------------------------------------------------------------- 
     
    189190      !---------------------------------------------------------------------- 
    190191 
    191                                                        ! =============== 
    192       DO jk = 1, jpkm1                                 ! Horizontal slab 
    193          !                                             ! =============== 
    194          ! Initialization of metric arrays (for z- or s-coordinates) 
    195          ! --------------------------------------------------------- 
    196          DO jj = 1, jpjm1 
    197             DO ji = 1, fs_jpim1   ! vector opt. 
    198 #if defined key_zco 
    199                ! z-coordinates, no vertical scale factors 
    200                zee(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk) 
    201 #else 
    202                ! vertical scale factor are used 
    203                zee(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 
    204 #endif 
    205             END DO 
    206          END DO 
     192      DO jk = 1, jpkm1 
    207193 
    208194         ! Laplacian of tracers (at before time step) 
     
    211197         DO jj = 1, jpjm1 
    212198            DO ji = 1, fs_jpim1   ! vector opt. 
    213                zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) ) 
    214             END DO 
    215          END DO 
    216          DO jj = 2, jpjm1 
    217             DO ji = fs_2, fs_jpim1   ! vector opt. 
    218 #if defined key_zco 
    219                zee(ji,jj) = e1t(ji,jj) /  e2t(ji,jj) 
    220 #else 
    221                zee(ji,jj) = e1t(ji,jj) / (e2t(ji,jj) * fse3t(ji,jj,jk)) 
    222 #endif 
    223                zlap(ji,jj,jk) = zee(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji-1,jj,jk)  ) 
     199               zmask(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk)                    & 
     200                  &                         * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) ) / e1u(ji,jj) 
     201            END DO 
     202         END DO 
     203         DO jj = 2, jpjm1 
     204            DO ji = fs_2, fs_jpim1   ! vector opt. 
     205               zlap(ji,jj,jk) = e1t(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji-1,jj,jk)  )   & 
     206                  &                        / ( e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    224207            END DO 
    225208         END DO 
     
    231214            DO ji = 2, fs_jpim1   ! vector opt. 
    232215               ! Upstream in the x-direction for the tracer 
    233                zmask(ji,jj,jk) =ptb(ji-1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji-1,jj,jk)) 
     216               zmask(ji,jj,jk) = ptb(ji-1,jj,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji-1,jj,jk) ) 
    234217               ! Downstream in the x-direction for the tracer 
    235                dwst (ji,jj,jk) =ptb(ji+1,jj,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji+1,jj,jk)) 
    236             ENDDO 
    237          ENDDO 
    238       END DO 
    239       !--- Lateral boundary conditions on the laplacian (unchanged sgn) 
    240       CALL lbc_lnk(   zlap(:,:,:), 'T', 1. )  
    241       !--- Lateral boundary conditions for the lim function 
    242       CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  dwst(:,:,:), 'T', 1. ) 
    243                                                        ! =============== 
    244       DO jk = 1, jpkm1                                 ! Horizontal slab 
    245          !                                             ! =============== 
    246          ! --- lim at the U-points in function of the direction of the flow 
    247          ! ---------------------------------------------------------------- 
     218               zdwst(ji,jj,jk) = ptb(ji+1,jj,jk) + sl(ji,jj,jk) * ( ptb(ji,jj,jk) - ptb(ji+1,jj,jk) ) 
     219            END DO 
     220         END DO 
     221      END DO 
     222 
     223      !--- Lateral boundary conditions on the laplacian and lim functions (unchanged sgn) 
     224      CALL lbc_lnk( zlap (:,:,:), 'T', 1. )  
     225      CALL lbc_lnk( zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  zdwst(:,:,:), 'T', 1. ) 
     226             
     227      ! --- lim at the U-points in function of the direction of the flow 
     228      ! ---------------------------------------------------------------- 
     229      DO jk = 1, jpkm1 
    248230         DO jj = 1, jpjm1 
    249231            DO ji = 1, fs_jpim1   ! vector opt. 
    250                dir = 0.5 + sign(0.5,pun(ji,jj,jk))     ! if pun>0 : diru = 1 otherwise diru = 0 
    251                lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji+1,jj,jk) 
     232               zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )     ! if pun>0 : diru = 1 otherwise diru = 0 
     233               zlim(ji,jj,jk) = zdir * zmask(ji,jj,jk) + (1-zdir) * zdwst(ji+1,jj,jk) 
    252234               ! Mask at the T-points in the x-direction (mask=0 or mask=1) 
    253                zmask(ji,jj,jk)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2 
    254             END DO 
    255          END DO 
    256       END DO 
    257       !--- Lateral boundary conditions for the mask 
    258       CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )   
     235               zmask(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2 
     236            END DO 
     237         END DO 
     238      END DO 
     239      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )       !--- Lateral boundary conditions for the mask 
     240 
    259241 
    260242      ! Horizontal advective fluxes 
    261243      ! --------------------------- 
    262                                                        ! =============== 
    263       DO jk = 1, jpkm1                                 ! Horizontal slab 
    264                                                        ! =============== 
    265          dt  = z2 * rdttra(jk) 
     244      DO jk = 1, jpkm1 
     245         zdt  = pz2 * rdttra(jk) 
    266246         !--- tracer flux at u and v-points 
    267247         DO jj = 1, jpjm1 
    268248            DO ji = 1, fs_jpim1   ! vector opt.  
    269 #if defined key_zco 
    270                e2 = e2u(ji,jj) 
    271 #else 
    272                e2 = e2u(ji,jj) * fse3u(ji,jj,jk) 
    273 #endif 
    274                dir = 0.5 + sign(0.5,pun(ji,jj,jk))                ! if pun>0 : diru = 1 otherwise diru = 0 
    275  
    276                dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj) 
    277                c  = ABS(pun(ji,jj,jk))*dt/dx                      ! (0<cx<1 : Courant number on x-direction) 
    278  
    279                fu  = lim(ji,jj,jk)                                ! FU + sl(FC-FU) in the x-direction for T 
    280                fc  = dir*ptb(ji  ,jj,jk)+(1-dir)*ptb(ji+1,jj,jk)  ! FC in the x-direction for T 
    281                fd  = dir*ptb(ji+1,jj,jk)+(1-dir)*ptb(ji  ,jj,jk)  ! FD in the x-direction for T 
     249               zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) )                     ! if pun>0 : diru = 1 otherwise diru = 0 
     250               zdx  = ( zdir * e1t(ji,jj) + (1-zdir)* e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk) 
     251               zc   = ABS( pun(ji,jj,jk) ) * zdt / zdx                     ! (0<cx<1 : Courant number on x-direction) 
     252 
     253               zfu  = zlim(ji,jj,jk)                                       ! FU + sl(FC-FU) in the x-direction for T 
     254               zfc  = zdir * ptb(ji  ,jj,jk) + (1-zdir) * ptb(ji+1,jj,jk)  ! FC in the x-direction for T 
     255               zfd  = zdir * ptb(ji+1,jj,jk) + (1-zdir) * ptb(ji  ,jj,jk)  ! FD in the x-direction for T 
    282256 
    283257               !--- QUICKEST scheme 
    284258               ! Temperature on the x-direction 
    285                coef1 = 0.5*(fc+fd) 
    286                coef2 = 0.5*c*(fd-fc) 
    287                coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji+1,jj,jk) ) 
    288                fho = coef1-coef2-coef3 
    289                fho = bound(fu,fd,fc,fho) 
     259               zcoef1 = 0.5      * ( zfc + zfd ) 
     260               zcoef2 = 0.5 * zc * ( zfd - zfc ) 
     261               zcoef3 = ( (1.-(zc*zc)) / 6. ) * ( zdir * zlap(ji,jj,jk) + (1-zdir) * zlap(ji+1,jj,jk) ) 
     262               zfho = zcoef1 - zcoef2 - zcoef3 
     263               zfho = bound( zfu, zfd, zfc, zfho ) 
    290264               !--- If the second ustream point is a land point 
    291265               !--- the flux is computed by the 1st order UPWIND scheme 
    292                mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji+1,jj,jk) 
    293                fho = mask*fho + (1-mask)*fc 
    294                dwst(ji,jj,jk)=e2*pun(ji,jj,jk)*fho 
     266               zmsk = zdir * zmask(ji,jj,jk) + (1-zdir) * zmask(ji+1,jj,jk) 
     267               zfho = zmsk * zfho            + (1-zmsk) * zfc 
     268               zdwst(ji,jj,jk) = pun(ji,jj,jk) * zfho 
    295269            END DO 
    296270         END DO 
     
    300274            DO ji = fs_2, fs_jpim1   ! vector opt. 
    301275               !--- horizontal advective trends 
    302 #if defined key_zco 
    303                zbtr = zbtr2(ji,jj) 
    304 #else 
    305                zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)               
    306 #endif 
    307                za = - zbtr * ( dwst(ji,jj,jk) - dwst(ji-1,jj  ,jk) ) 
     276               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    308277               !--- add it to the general tracer trends 
    309                pta(ji,jj,jk) = pta(ji,jj,jk) + za 
    310             END DO 
    311          END DO 
    312          !                                             ! =============== 
    313       END DO                                           !   End of slab 
    314       !                                                ! =============== 
    315  
    316       ! Save the horizontal advective trends for diagnostic 
    317       IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, dwst, pun, ptb ) 
     278               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * ( zdwst(ji,jj,jk) - zdwst(ji-1,jj  ,jk) ) 
     279            END DO 
     280         END DO 
     281      END DO 
     282 
     283      !                                            !-- trend diagnostic 
     284      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zdwst, pun, ptb ) 
    318285 
    319286 
     
    321288      ! I. Part 2 : y-direction 
    322289      !---------------------------------------------------------------------- 
    323                                                        ! ============== 
    324       DO jk = 1, jpkm1                                 ! Horizontal slab 
    325          !                                             ! =============== 
    326          ! Initialization of metric arrays (for z- or s-coordinates) 
    327          ! --------------------------------------------------------- 
    328          DO jj = 1, jpjm1 
    329             DO ji = 1, fs_jpim1   ! vector opt. 
    330 #if defined key_zco 
    331                ! z-coordinates, no vertical scale factors 
    332                zee(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk) 
    333 #else 
    334                ! s-coordinates, vertical scale factor are used 
    335                zee(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 
    336 #endif 
    337             END DO 
    338          END DO 
     290 
     291      DO jk = 1, jpkm1 
    339292 
    340293         ! Laplacian of tracers (at before time step) 
     
    343296         DO jj = 1, jpjm1 
    344297            DO ji = 1, fs_jpim1   ! vector opt. 
    345                zmask(ji,jj,jk) = zee(ji,jj) * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) ) 
     298               zmask(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk)   & 
     299                  &            * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk) 
    346300            END DO 
    347301         END DO 
     
    349303         DO jj = 2, jpjm1 
    350304            DO ji = fs_2, fs_jpim1   ! vector opt. 
    351 #if defined key_zco 
    352                zee(ji,jj) = e2t(ji,jj) /  e1t(ji,jj) 
    353 #else 
    354                zee(ji,jj) = e2t(ji,jj) / (e1t(ji,jj) * fse3t(ji,jj,jk)) 
    355 #endif 
    356                zlap(ji,jj,jk) = zee(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji,jj-1,jk)  ) 
     305               zlap(ji,jj,jk) = e2t(ji,jj) * (  zmask(ji,jj,jk) - zmask(ji,jj-1,jk)  )   & 
     306                  &                        / ( e1t(ji,jj) * fse3t(ji,jj,jk) ) 
    357307            END DO 
    358308         END DO 
     
    362312            DO ji = 2, fs_jpim1   ! vector opt. 
    363313               ! Upstream in the y-direction for the tracer 
    364                zmask(ji,jj,jk)=ptb(ji,jj-1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj-1,jk)) 
     314               zmask(ji,jj,jk) = ptb(ji,jj-1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj-1,jk) ) 
    365315               ! Downstream in the y-direction for the tracer 
    366                dwst (ji,jj,jk)=ptb(ji,jj+1,jk)+sl(ji,jj,jk)*(ptb(ji,jj,jk)-ptb(ji,jj+1,jk)) 
     316               zdwst(ji,jj,jk) = ptb(ji,jj+1,jk) + sl(ji,jj,jk) *( ptb(ji,jj,jk) - ptb(ji,jj+1,jk) ) 
    367317            ENDDO 
    368318         ENDDO 
    369319      END DO 
    370       !--- Lateral boundary conditions on the laplacian (unchanged sgn) 
    371       CALL lbc_lnk(   zlap(:,:,:), 'T', 1. ) 
    372       !--- Lateral boundary conditions for the lim function 
    373       CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk(  dwst(:,:,:), 'T', 1. ) 
    374  
    375       DO jk = 1, jpkm1                                 ! Horizontal slab 
    376          !                                             ! =============== 
     320      !--- Lateral boundary conditions on the laplacian and lim function (unchanged sgn) 
     321      CALL lbc_lnk( zlap (:,:,:), 'T', 1. ) 
     322      CALL lbc_lnk( zmask(:,:,:), 'T', 1. )   ;   CALL lbc_lnk( zdwst(:,:,:), 'T', 1. ) 
     323 
     324      DO jk = 1, jpkm1 
     325       
    377326         ! --- lim at the V-points in function of the direction of the flow 
    378327         ! ---------------------------------------------------------------- 
    379328         DO jj = 1, jpjm1 
    380329            DO ji = 1, fs_jpim1   ! vector opt. 
    381                dir = 0.5 + sign(0.5,pvn(ji,jj,jk))      ! if pvn>0 : dirv = 1 otherwise dirv = 0 
    382                lim(ji,jj,jk)=dir*zmask(ji,jj,jk)+(1-dir)*dwst(ji,jj+1,jk) 
     330               zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )      ! if pvn>0 : dirv = 1 otherwise dirv = 0 
     331               zlim(ji,jj,jk) = zdir * zmask(ji,jj,jk) + (1-zdir) * zdwst(ji,jj+1,jk) 
    383332               ! Mask at the T-points in the y-direction (mask=0 or mask=1) 
    384                zmask(ji,jj,jk)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2 
    385             END DO 
    386          END DO 
    387       END DO 
    388       !--- Lateral boundary conditions for the mask 
    389       CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )  
     333               zmask(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 
     334            END DO 
     335         END DO 
     336      END DO 
     337      CALL lbc_lnk(  zmask(:,:,:), 'T', 1. )      !--- Lateral boundary conditions for the mask 
    390338 
    391339      ! Horizontal advective fluxes 
    392       ! ------------------------------- 
    393                                                        ! =============== 
    394       DO jk = 1, jpkm1                                 ! Horizontal slab 
    395                                                        ! =============== 
    396          dt  = z2 * rdttra(jk) 
     340      ! --------------------------- 
     341                                   
     342      DO jk = 1, jpkm1            
     343         zdt  = pz2 * rdttra(jk) 
    397344         !--- tracer flux at u and v-points 
    398345         DO jj = 1, jpjm1 
    399346            DO ji = 1, fs_jpim1   ! vector opt. 
    400 #if defined key_zco 
    401                e1 = e1v(ji,jj) 
    402 #else 
    403                e1 = e1v(ji,jj) * fse3v(ji,jj,jk) 
    404 #endif 
    405                dir = 0.5 + sign(0.5,pvn(ji,jj,jk))                ! if pvn>0 : dirv = 1 otherwise dirv = 0 
    406  
    407                dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1) 
    408                c  = ABS(pvn(ji,jj,jk))*dt/dx                      ! (0<cy<1 : Courant number on y-direction) 
    409  
    410                fu  = lim(ji,jj,jk)                                ! FU + sl(FC-FU) in the y-direction for T 
    411                fc  = dir*ptb(ji,jj  ,jk)+(1-dir)*ptb(ji,jj+1,jk)  ! FC in the y-direction for T 
    412                fd  = dir*ptb(ji,jj+1,jk)+(1-dir)*ptb(ji,jj  ,jk)  ! FD in the y-direction for T 
     347               zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) )                     ! if pvn>0 : dirv = 1 otherwise dirv = 0 
     348               zdx  = ( zdir * e2t(ji,jj) + (1-zdir)* e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk) 
     349               zc   = ABS( pvn(ji,jj,jk) ) * zdt / zdx                     ! (0<cy<1 : Courant number on y-direction) 
     350 
     351               zfu  = zlim(ji,jj,jk)                                       ! FU + sl(FC-FU) in the y-direction for T 
     352               zfc  = zdir * ptb(ji,jj  ,jk) + (1-zdir) * ptb(ji,jj+1,jk)  ! FC in the y-direction for T 
     353               zfd  = zdir * ptb(ji,jj+1,jk) + (1-zdir) * ptb(ji,jj  ,jk)  ! FD in the y-direction for T 
    413354 
    414355               !--- QUICKEST scheme 
    415356               ! Temperature on the y-direction 
    416                coef1 = 0.5*(fc+fd) 
    417                coef2 = 0.5*c*(fd-fc) 
    418                coef3 = ((1.-(c*c))/6.)*(dir*zlap(ji,jj,jk) + (1-dir)*zlap(ji,jj+1,jk) ) 
    419                fho = coef1-coef2-coef3 
    420                fho = bound(fu,fd,fc,fho) 
     357               zcoef1 = 0.5      * ( zfc + zfd ) 
     358               zcoef2 = 0.5 * zc * ( zfd - zfc ) 
     359               zcoef3 = ( (1.-(zc*zc)) / 6. ) * ( zdir * zlap(ji,jj,jk) + (1-zdir) * zlap(ji,jj+1,jk) ) 
     360               zfho = zcoef1 - zcoef2 - zcoef3 
     361               zfho = bound( zfu, zfd, zfc, zfho ) 
    421362               !--- If the second ustream point is a land point 
    422363               !--- the flux is computed by the 1st order UPWIND scheme 
    423                mask=dir*zmask(ji,jj,jk)+(1-dir)*zmask(ji,jj+1,jk) 
    424                fho = mask*fho + (1-mask)*fc 
    425                dwst(ji,jj,jk)=e1*pvn(ji,jj,jk)*fho 
     364               zmsk = zdir * zmask(ji,jj,jk) + (1-zdir) * zmask(ji,jj+1,jk) 
     365               zfho = zmsk * zfho            + (1-zmsk) * zfc 
     366               zdwst(ji,jj,jk) = pvn(ji,jj,jk) * zfho 
    426367            END DO 
    427368         END DO 
     
    430371         DO jj = 2, jpjm1 
    431372            DO ji = fs_2, fs_jpim1   ! vector opt. 
    432                !--- horizontal advective trends 
    433 #if defined key_zco 
    434                zbtr = zbtr2(ji,jj) 
    435 #else 
    436                zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk) 
    437 #endif 
    438                za = - zbtr * (  dwst(ji,jj,jk) - dwst(ji  ,jj-1,jk) ) 
    439                !--- add it to the general tracer trends 
    440                pta(ji,jj,jk) = pta(ji,jj,jk) + za 
    441             END DO 
    442          END DO 
    443          !                                             ! =============== 
    444       END DO                                           !   End of slab 
    445       !                                                ! =============== 
    446  
    447       ! Save the horizontal advective trends for diagnostic 
    448       IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, dwst, pvn, ptb ) 
    449  
    450       ! "zonal" mean advective heat and salt transport  
     373               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     374               ! horizontal advective trends added to the general tracer trends 
     375               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * (  zdwst(ji,jj,jk) - zdwst(ji,jj-1,jk) ) 
     376            END DO 
     377         END DO 
     378      END DO  
     379 
     380      !                                           !-- trend diagnostic 
     381      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zdwst, pvn, ptb ) 
     382      !                                           ! "Poleward" heat or salt transport  
    451383      IF( ln_diaptr .AND. cdtype == 'TRA' .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
    452 #if defined key_zco 
    453          DO jk = 1, jpkm1 
    454             DO jj = 2, jpjm1 
    455                DO ji = fs_2, fs_jpim1   ! vector opt. 
    456                   dwst(ji,jj,jk) = dwst(ji,jj,jk) * fse3v(ji,jj,jk) 
    457                END DO 
    458             END DO 
    459          END DO 
    460 # endif 
    461          IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( dwst(:,:,:) ) 
    462          IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( dwst(:,:,:) ) 
     384         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zdwst(:,:,:) ) 
     385         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zdwst(:,:,:) ) 
    463386      ENDIF 
    464387      ! 
     
    466389 
    467390 
    468    SUBROUTINE tra_adv_qck_ver( pwn, ptn , pta, z2  )       
    469       !!---------------------------------------------------------------------- 
    470       !! 
    471       !!---------------------------------------------------------------------- 
    472       REAL(wp), INTENT(in   )                         ::   z2 
    473       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn 
    474       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptn 
    475       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta 
    476       !! 
    477       REAL(wp) ::   za, ze3tr, dt, dir, fc, fd             ! temporary scalars 
     391   SUBROUTINE tra_adv_qck_ver( pwn, ptn , pta )       
     392      !!---------------------------------------------------------------------- 
     393      !!                  ***  ROUTINE tra_adv_qck_ver  *** 
     394      !! 
     395      !! ** Purpose :    
     396      !! 
     397      !!---------------------------------------------------------------------- 
     398      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn   ! vertical transport 
     399      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptn   ! now tracer 
     400      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta   ! tracer trend 
     401      !! 
     402      INTEGER  ::   ji, jj , jk                  ! dummy loop indices 
     403      REAL(wp) ::   zbtr, zdir, zfc, zfd   ! temporary scalars 
     404      !!---------------------------------------------------------------------- 
    478405 
    479406      ! Vertical advection 
     
    483410      ! ---------------------------- 
    484411 
    485       sl(:,:,jpk) = 0.e0      !Bottom value : flux set to zero 
    486  
    487       ! Surface value 
    488       IF( lk_dynspg_rl .OR. lk_vvl ) THEN         ! rigid lid : flux set to zero 
    489          sl(:,:, 1 ) = 0.e0 
    490       ELSE                                        ! free surface-constant volume 
    491          sl(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) 
     412      sl(:,:,jpk) = 0.e0                          ! bottom value  
     413 
     414      !                                           ! Surface value 
     415      IF( lk_dynspg_rl .OR. lk_vvl ) THEN   ;   sl(:,:, 1 ) = 0.e0                      ! rigid lid or non-linear fs 
     416      ELSE                                  ;   sl(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1)   ! linear free surface 
    492417      ENDIF 
     418 
     419!!gm  bug: code au true 2nd order scheme 
     420!!gm  : sl used as work array: not good idea for optimisation (sl compute once for all tracers...) 
    493421 
    494422      ! Second order centered tracer flux at w-point 
    495423      DO jk = 2, jpkm1 
    496          dt  = z2 *  rdttra(jk) 
    497          DO jj = 2, jpjm1 
    498             DO ji = fs_2, fs_jpim1   ! vector opt. 
    499                dir = 0.5 + sign(0.5,pwn(ji,jj,jk))                                         ! if pwn>0 : dirw = 1 otherwise dirw = 0 
    500                fc = dir*ptn(ji,jj,jk  )*fse3t(ji,jj,jk-1)+(1-dir)*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk  )   ! FC in the z-direction for T 
    501                fd = dir*ptn(ji,jj,jk-1)*fse3t(ji,jj,jk  )+(1-dir)*ptn(ji,jj,jk  )*fse3t(ji,jj,jk-1)   ! FD in the z-direction for T 
     424         DO jj = 2, jpjm1 
     425            DO ji = fs_2, fs_jpim1   ! vector opt. 
     426               zdir = 0.5 + SIGN( 0.5, pwn(ji,jj,jk) )                        ! if pwn>0 : dirw = 1 otherwise dirw = 0 
     427               zfc = zdir * ptn(ji,jj,jk  ) * fse3t(ji,jj,jk-1) + (1-zdir) * ptn(ji,jj,jk-1) * fse3t(ji,jj,jk  )   ! FC  
     428               zfd = zdir * ptn(ji,jj,jk-1) * fse3t(ji,jj,jk  ) + (1-zdir) * ptn(ji,jj,jk  ) * fse3t(ji,jj,jk-1)   ! FD 
    502429               !--- Second order centered scheme 
    503                sl(ji,jj,jk)=pwn(ji,jj,jk)*(fc+fd)/(fse3t(ji,jj,jk-1)+fse3t(ji,jj,jk)) 
     430               sl(ji,jj,jk) = pwn(ji,jj,jk) * ( zfc + zfd ) / ( fse3t(ji,jj,jk-1) + fse3t(ji,jj,jk) ) 
    504431            END DO 
    505432         END DO 
     
    511438         DO jj = 2, jpjm1 
    512439            DO ji = fs_2, fs_jpim1   ! vector opt. 
    513                ze3tr = 1. / fse3t(ji,jj,jk) 
    514                ! vertical advective trends 
    515                za = - ze3tr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) ) 
    516                ! add it to the general tracer trends 
    517                pta(ji,jj,jk) =  pta(ji,jj,jk) + za 
     440               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     441               ! vertical advective trends added to the general tracer trends 
     442               pta(ji,jj,jk) =  pta(ji,jj,jk) - zbtr * ( sl(ji,jj,jk) - sl(ji,jj,jk+1) ) 
    518443            END DO 
    519444         END DO 
     
    523448 
    524449 
    525    REAL FUNCTION bound(fu,fd,fc,fho) 
    526       REAL(wp) ::  fu, fd, fc, fho, fref1, fref2 
    527       fref1 = fu 
    528       fref2 = MAX(  MIN( fc , fd ), MIN( MAX( fc , fd ), fref1 )  ) 
    529       bound = MAX(  MIN( fho, fc ), MIN( MAX( fho, fc ), fref2 )  ) 
     450   FUNCTION bound( pfu, pfd, pfc, pfho )   RESULT( pbound ) 
     451      !!---------------------------------------------------------------------- 
     452      !!                  ***  FUNCTION bound  *** 
     453      !! 
     454      !! ** Purpose :   ??? 
     455      !!---------------------------------------------------------------------- 
     456      REAL(wp), INTENT(in) ::  pfu, pfd, pfc, pfho 
     457      REAL(wp)             ::  pbound 
     458      !! 
     459      REAL(wp) ::  zfref1, zfref2 
     460      !!---------------------------------------------------------------------- 
     461      zfref1 = pfu 
     462      zfref2 = MAX(  MIN( pfc , pfd ), MIN( MAX( pfc , pfd ), zfref1 )  ) 
     463      pbound = MAX(  MIN( pfho, pfc ), MIN( MAX( pfho, pfc ), zfref2 )  ) 
     464      ! 
    530465   END FUNCTION 
    531466 
  • branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r786 r791  
    44   !! Ocean active tracers:  horizontal & vertical advective trend 
    55   !!============================================================================== 
    6    !! History :  7.0  !  95-12  (L. Mortier)  Original code 
    7    !!            8.0  !  00-01  (H. Loukos)  adapted to ORCA  
    8    !!             -   !  00-10  (MA Foujols E.Kestenare)  include file not routine 
    9    !!             -   !  00-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes 
    10    !!             -   !  01-07  (E. Durand G. Madec)  adaptation to ORCA config 
    11    !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
    12    !!   NEMO     1.0  !  04-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl 
    13    !!             -   !  08-04  (S. Cravatte) add the i-, j- & k- trends computation 
    14    !!             -   !  05-11  (V. Garnier) Surface pressure gradient organization 
    15    !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC 
     6   !! History :  7.0  !  1995-12  (L. Mortier)  Original code 
     7   !!            8.0  !  2000-01  (H. Loukos)  adapted to ORCA  
     8   !!             -   !  2000-10  (MA Foujols E.Kestenare)  include file not routine 
     9   !!             -   !  2000-12  (E. Kestenare M. Levy)  fix bug in trtrd indexes 
     10   !!             -   !  2001-07  (E. Durand G. Madec)  adaptation to ORCA config 
     11   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module 
     12   !!   NEMO     1.0  !  2004-01  (A. de Miranda, G. Madec, J.M. Molines ): advective bbl 
     13   !!             -   !  2008-04  (S. Cravatte) add the i-, j- & k- trends computation 
     14   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization 
     15   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport 
    1616   !!---------------------------------------------------------------------- 
    1717 
     
    4242   !!---------------------------------------------------------------------- 
    4343   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)  
    44    !! $Id:$  
     44   !! $Id$  
    4545   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4646   !!---------------------------------------------------------------------- 
     
    7070      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend  
    7171      !! 
    72       INTEGER  ::   ji, jj, jk              ! dummy loop indices 
    73       REAL(wp) ::   ztai, ztaj, ztak 
    74       REAL(wp) ::   z2dtt, zbtr, zeu, zev   ! temporary scalar 
    75       REAL(wp) ::   zew, z2                          ! temporary scalar 
    76       REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk           !    "         " 
    77       REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk            !    "         " 
    78       REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztu, ztv, ztw   ! temporary workspace 
     72      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
     73      REAL(wp) ::   z2dtt, zbtr, z2, zzti    ! temporary scalar 
     74      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !    "         " 
     75      REAL(wp) ::   zfm_ui, zfm_vj, zfm_wk   !    "         " 
     76      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztu, ztv, ztw   ! 3D workspace 
    7977      !!---------------------------------------------------------------------- 
    8078 
     
    8785      ENDIF 
    8886 
    89       IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1. 
    90       ELSE                                        ;    z2 = 2. 
     87      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2 = 1.      ! euler   time-stepping 
     88      ELSE                                        ;    z2 = 2.      ! leap-frog time-stepping 
    9189      ENDIF 
    9290 
     
    103101         DO jj = 1, jpjm1 
    104102            DO ji = 1, fs_jpim1   ! vector opt. 
    105                zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 
    106                zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 
    107                ! upstream scheme 
    108                zfp_ui = zeu + ABS( zeu ) 
    109                zfm_ui = zeu - ABS( zeu ) 
    110                zfp_vj = zev + ABS( zev ) 
    111                zfm_vj = zev - ABS( zev ) 
    112                ztu(ji,jj,jk) = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk) 
    113                ztv(ji,jj,jk) = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk) 
    114             END DO 
    115          END DO 
    116       END DO 
    117  
     103               zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 
     104               zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 
     105               zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 
     106               zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 
     107               ztu(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk) ) 
     108               ztv(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk) ) 
     109            END DO 
     110         END DO 
     111      END DO 
     112      ! 
    118113      ! upstream tracer flux in the k direction 
    119       ! Surface value 
    120       IF( lk_dynspg_rl .OR. lk_vvl ) THEN               ! rigid lid or variable volume: flux set to zero 
    121          ztw(:,:,1) = 0.e0 
    122       ELSE                                              ! free surface 
    123          ztw(:,:,1) = e1t(:,:) * e2t(:,:) * pwn(:,:,1) * ptb(:,:,1) 
    124       ENDIF 
    125  
    126       ! Interior value 
    127       DO jk = 2, jpkm1 
     114      !                                                                           ! Surface value 
     115      IF( lk_dynspg_rl .OR. lk_vvl ) THEN   ;   ztw(:,:,1) = 0.e0                      ! rigid lid or non-linear fs 
     116      ELSE                                  ;   ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1)   ! free surface 
     117      ENDIF 
     118      ! 
     119      DO jk = 2, jpkm1                                                            ! Interior value 
    128120         DO jj = 1, jpj 
    129121            DO ji = 1, jpi 
    130                zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 
    131                zfp_wk = zew + ABS( zew ) 
    132                zfm_wk = zew - ABS( zew ) 
    133                ztw(ji,jj,jk) = zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1) 
    134             END DO 
    135          END DO 
    136       END DO 
    137  
    138       ! total advective trend 
    139       DO jk = 1, jpkm1 
     122               zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     123               zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     124               ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1) ) 
     125            END DO 
     126         END DO 
     127      END DO 
     128      ! 
     129      ! upstream advective trend 
     130      DO jk = 1, jpkm1 
     131         z2dtt = z2 * rdttra(jk) 
    140132         DO jj = 2, jpjm1 
    141133            DO ji = fs_2, fs_jpim1   ! vector opt. 
    142134               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    143                ! i- j- horizontal & k- vertical advective trends 
    144                ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  ) ) * zbtr 
    145                ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  ) ) * zbtr 
    146                ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr 
    147135               ! total intermediate advective trends 
    148                zti(ji,jj,jk) = ztai + ztaj + ztak 
    149             END DO 
    150          END DO 
    151       END DO 
    152  
    153       !  Save the horizontal advective trends for diagnostic 
    154       ! ----------------------------------------------------- 
     136               zzti = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )            & 
     137                  &     + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )            & 
     138                  &     + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) * zbtr 
     139               ! update and guess with monotonic sheme 
     140               pta(ji,jj,jk) =   pta(ji,jj,jk) +         zzti 
     141               zti(ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * zzti ) * tmask(ji,jj,jk) 
     142            END DO 
     143         END DO 
     144      END DO 
     145      ! 
     146      CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti   (unchanged sign) 
     147 
     148 
     149      !                                 ! trend diagnostics (contribution of upstream fluxes) 
    155150      IF( l_trdtra ) THEN 
    156151         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn ) 
     
    158153         CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn ) 
    159154      ENDIF 
    160  
    161       ! update and guess with monotonic sheme 
    162       DO jk = 1, jpkm1 
    163          z2dtt = z2 * rdttra(jk) 
    164          DO jj = 2, jpjm1 
    165             DO ji = fs_2, fs_jpim1   ! vector opt. 
    166                pta(ji,jj,jk) =  pta(ji,jj,jk) + zti(ji,jj,jk) 
    167                zti (ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * zti(ji,jj,jk) ) * tmask(ji,jj,jk) 
    168             END DO 
    169          END DO 
    170       END DO 
    171  
    172       ! Lateral boundary conditions on zti   (unchanged sign) 
    173       CALL lbc_lnk( zti, 'T', 1. ) 
     155      !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     156      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
     157         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( ztv(:,:,:) )  
     158         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( ztv(:,:,:) ) 
     159      ENDIF 
    174160 
    175161 
    176162      ! 3. antidiffusive flux : high order minus low order 
    177163      ! -------------------------------------------------- 
    178       ! antidiffusive flux on i and j 
     164      !                                 ! anti-diffusive flux on i and j 
    179165      DO jk = 1, jpkm1 
    180166         DO jj = 1, jpjm1 
    181167            DO ji = 1, fs_jpim1   ! vector opt. 
    182                zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 
    183                zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 
    184                ztu(ji,jj,jk) = zeu * ( ptn(ji,jj,jk) + ptn(ji+1,jj,jk) ) - ztu(ji,jj,jk) 
    185                ztv(ji,jj,jk) = zev * ( ptn(ji,jj,jk) + ptn(ji,jj+1,jk) ) - ztv(ji,jj,jk) 
    186             END DO 
    187          END DO 
    188       END DO 
    189        
    190       ! antidiffusive flux on k 
    191       ! Surface value 
    192       ztw(:,:,1) = 0.e0 
    193  
    194       ! Interior value 
    195       DO jk = 2, jpkm1 
     168               ztu(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj,jk) ) - ztu(ji,jj,jk) 
     169               ztv(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj+1,jk) ) - ztv(ji,jj,jk) 
     170            END DO 
     171         END DO 
     172      END DO 
     173      !                                 ! antidiffusive flux on k 
     174      ztw(:,:,1) = 0.e0                      ! Surface value 
     175      ! 
     176      DO jk = 2, jpkm1                       ! Interior value 
    196177         DO jj = 1, jpj 
    197178            DO ji = 1, jpi 
    198                zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 
    199                ztw(ji,jj,jk) = zew * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 
    200             END DO 
    201          END DO 
    202       END DO 
    203  
    204       ! Lateral bondary conditions 
    205       CALL lbc_lnk( ztu, 'U', -1. ) 
     179               ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 
     180            END DO 
     181         END DO 
     182      END DO 
     183      ! 
     184      CALL lbc_lnk( ztu, 'U', -1. )     ! Lateral bondary conditions on upstream fluxes 
    206185      CALL lbc_lnk( ztv, 'V', -1. ) 
    207186      CALL lbc_lnk( ztw, 'W',  1. ) 
    208187 
    209       ! 4. monotonicity algorithm 
    210       ! ------------------------- 
     188      !                                 ! monotonicity algorithm 
    211189      CALL nonosc( ptb, ztu, ztv, ztw, zti, z2 ) 
    212190 
    213191 
    214       ! 5. final trend with corrected fluxes 
    215       ! ------------------------------------ 
     192      ! 4. final trend with anti-diffusive fluxes 
     193      ! ----------------------------------------- 
    216194      DO jk = 1, jpkm1 
    217195         DO jj = 2, jpjm1 
    218196            DO ji = fs_2, fs_jpim1   ! vector opt.   
    219197               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    220                ! i- j- horizontal & k- vertical advective trends 
    221                ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )) * zbtr 
    222                ztaj = - ( ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )) * zbtr 
    223                ztak = - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1)) * zbtr 
    224  
    225                ! add them to the general tracer trends 
    226                pta(ji,jj,jk) = pta(ji,jj,jk) + ztai + ztaj + ztak 
    227             END DO 
    228          END DO 
    229       END DO 
    230  
    231 !!gm  the transport computation is wrong, the upstream part is missing ! 
    232       ! "zonal" mean advective heat and salt transport 
    233       IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
    234          IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( ztv(:,:,:) ) 
    235          IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( ztv(:,:,:) ) 
    236       ENDIF 
    237  
    238       !  Save the horizontal advective trends for diagnostic 
    239       ! ----------------------------------------------------- 
    240       IF( l_trdtra ) THEN 
     198               ! anti-diffusive trends added to the general tracer trends 
     199               pta(ji,jj,jk) = pta(ji,jj,jk) - ( ztu(ji,jj,jk) - ztu(ji-1,jj  ,jk  )             & 
     200                  &                            + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk  )             & 
     201                  &                            + ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1)  ) * zbtr 
     202            END DO 
     203         END DO 
     204      END DO 
     205 
     206      IF( l_trdtra ) THEN               ! trend diagnostic (contribution of anti-diffusive fluxes) 
    241207         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, ztu, pun, ptn, cnbpas='bis' )   ! <<< Add to iad trend 
    242208         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, ztv, pvn, ptn, cnbpas='bis' )   ! <<< Add to jad trend 
    243209         CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, ztw, pwn, ptn, cnbpas='bis' )   ! <<< Add to zad trend 
    244210      ENDIF 
    245  
     211      !                                 ! "Poleward" heat and salt transports (contribution of anti-diffusive fluxes) 
     212      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
     213         IF( ktra == jp_tem)   pht_adv(:) = pht_adv(:) + ptr_vj( ztv(:,:,:) ) 
     214         IF( ktra == jp_sal)   pst_adv(:) = pst_adv(:) + ptr_vj( ztv(:,:,:) ) 
     215      ENDIF 
     216      !                                 ! control print 
    246217      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' tvd - adv: ', mask1=tmask, clinfo3=cdtype ) 
    247218      ! 
     
    276247      zbetup(:,:,:) = 0.e0   ;   zbetdo(:,:,:) = 0.e0   
    277248 
     249!!gm   optimisation :  add the optimal version I wrote 1 year ago 
    278250      ! Search local extrema 
    279251      ! -------------------- 
  • branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r786 r791  
    44   !! Ocean tracers:  horizontal & vertical advective trend 
    55   !!============================================================================== 
    6    !! History :  1.0  !  06-08  (L. Debreu, R. Benshila)  Original code 
    7    !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC 
     6   !! History :  1.0  !  2006-08  (L. Debreu, R. Benshila)  Original code 
     7   !!            2.4  !  2008-01  (G. Madec)  merge TRC-TRA + switch from velocity to transport 
    88   !!---------------------------------------------------------------------- 
    99 
     
    2727   PUBLIC   tra_adv_ubs   ! routine called by traadv module 
    2828 
    29    REAL(wp), DIMENSION(jpi,jpj) ::   e1e2tr   ! = 1/(e1t * e2t) 
    30  
    3129   !! * Substitutions 
    3230#  include "domzgr_substitute.h90" 
     
    3432   !!---------------------------------------------------------------------- 
    3533   !! NEMO/OPA & TRP 2.4 , LOCEAN-IPSL (2008)  
    36    !! $Id:$  
     34   !! $Id$  
    3735   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
    3836   !!---------------------------------------------------------------------- 
     
    7169      !! 
    7270      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.  
    73       !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 17311741.  
     71      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731-1741.  
    7472      !!---------------------------------------------------------------------- 
    7573      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index 
     
    8078      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend  
    8179      !! 
    82       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    83       REAL(wp) ::   zta, zbtr, zcoef                  ! temporary scalars 
    84       REAL(wp) ::   zfui, zfp_ui, zfm_ui, zcenut      !    "         " 
    85       REAL(wp) ::   zfvj, zfp_vj, zfm_vj, zcenvt      !    "         " 
    86       REAL(wp) ::   z_hdivn                           !    "         " 
    87       REAL(wp), DIMENSION(jpi,jpj)     ::   zeeu, zeev     ! temporary 2D workspace 
    88       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy              ! temporary 3D workspace 
    89       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu , ztv , zltu , zltv, ztrdt   !    "              " 
    90       !!---------------------------------------------------------------------- 
    91  
     80      INTEGER  ::   ji, jj, jk                     ! dummy loop indices 
     81      REAL(wp) ::   zbtr, zcoef, z_hdivn           ! temporary scalars 
     82      REAL(wp) ::   zeeu, zfp_ui, zfm_ui, zcenut   !    "         " 
     83      REAL(wp) ::   zeev, zfp_vj, zfm_vj, zcenvt   !    "         " 
     84      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, ztu, zltu   ! temporary 3D workspace 
     85      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwy, ztv, zltv   !    "              " 
     86      !!---------------------------------------------------------------------- 
     87 
     88!!gm  this can be optimized: we don't need zero everywhere 
    9289      zltu(:,:,:) = 0.e0 
    9390      zltv(:,:,:) = 0.e0 
     
    9794         IF(lwp) WRITE(numout,*) 'tra_adv_ubs :  horizontal UBS advection scheme' 
    9895         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    99          ! 
    100          e1e2tr(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
    101       ENDIF 
    102  
    103       ! store pta trends 
    104       ztrdt(:,:,:) = pta(:,:,:) 
    105  
    106       zcoef = 1./6. 
    107       !                                                ! =============== 
    108       DO jk = 1, jpkm1                                 ! Horizontal slab 
    109          !                                             ! =============== 
    110  
    111          !  Initialization of metric arrays (for z- or s-coordinates) 
     96      ENDIF 
     97 
     98      !  Laplacian 
     99      DO jk = 1, jpkm1 
     100         DO jj = 1, jpjm1                               ! First derivative (gradient) 
     101            DO ji = 1, fs_jpim1   ! vector opt. 
     102               zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 
     103               zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 
     104               ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) ) 
     105               ztv(ji,jj,jk) = zeev * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) ) 
     106            END DO 
     107         END DO 
     108         DO jj = 2, jpjm1                               ! Second derivative (divergence) 
     109            DO ji = fs_2, fs_jpim1   ! vector opt. 
     110               zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) 
     111               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef 
     112               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef 
     113            END DO 
     114         END DO 
     115      END DO 
     116 
     117      ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn) 
     118      CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. ) 
     119 
     120 
     121      !  Horizontal advective fluxes 
     122      DO jk = 1, jpkm1 
    112123         DO jj = 1, jpjm1 
    113124            DO ji = 1, fs_jpim1   ! vector opt. 
    114 #if defined key_zco 
    115                ! z-coordinates, no vertical scale factors 
    116                zeeu(ji,jj) = e2u(ji,jj) / e1u(ji,jj) * umask(ji,jj,jk) 
    117                zeev(ji,jj) = e1v(ji,jj) / e2v(ji,jj) * vmask(ji,jj,jk) 
    118 #else 
    119                ! s-coordinates, vertical scale factor are used 
    120                zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk) 
    121                zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk) 
    122 #endif 
    123             END DO 
    124          END DO 
    125  
    126          !  Laplacian 
    127          ! First derivative (gradient) 
    128          DO jj = 1, jpjm1 
    129             DO ji = 1, fs_jpim1   ! vector opt. 
    130                ztu(ji,jj,jk) = zeeu(ji,jj) * ( ptb(ji+1,jj  ,jk) - ptb(ji,jj,jk) ) 
    131                ztv(ji,jj,jk) = zeev(ji,jj) * ( ptb(ji  ,jj+1,jk) - ptb(ji,jj,jk) ) 
    132             END DO 
    133          END DO 
    134          ! Second derivative (divergence) 
    135          DO jj = 2, jpjm1 
    136             DO ji = fs_2, fs_jpim1   ! vector opt. 
    137 #if ! defined key_zco 
    138                zcoef = 1. / ( 6. * fse3t(ji,jj,jk) ) 
    139 #endif          
    140                zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef 
    141                zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef 
    142             END DO 
    143          END DO 
    144          !                                             ! ================= 
    145       END DO                                           !    End of slab 
    146       !                                                ! ================= 
    147  
    148       ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn) 
    149       CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. ) 
    150  
    151       !                                                ! =============== 
    152       DO jk = 1, jpkm1                                 ! Horizontal slab 
    153          !                                             ! =============== 
    154          !  Horizontal advective fluxes 
    155          DO jj = 1, jpjm1 
    156             DO ji = 1, fs_jpim1   ! vector opt. 
    157                ! volume fluxes * 1/2 
    158 #if defined key_zco 
    159                zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk) 
    160                zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk) 
    161 #else 
    162                zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk) 
    163                zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk) 
    164 #endif 
    165                ! upstream scheme 
    166                zfp_ui = zfui + ABS( zfui ) 
    167                zfp_vj = zfvj + ABS( zfvj ) 
    168                zfm_ui = zfui - ABS( zfui ) 
    169                zfm_vj = zfvj - ABS( zfvj ) 
     125               ! upstream transport 
     126               zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )      ! = 2.pun  if pun>0,  =0      if pun<0 
     127               zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )      ! = 0      if pun>0,  = 2.pun if pun<0 
     128               zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )      ! idem on pvn 
     129               zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )      ! 
    170130               ! centered scheme 
    171                zcenut = zfui * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) ) 
    172                zcenvt = zfvj * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) ) 
    173                ! mixed centered / upstream scheme 
    174                zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk) 
    175                zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk) 
    176             END DO 
    177          END DO 
    178  
    179          !  Tracer flux divergence at t-point added to the general trend 
    180          DO jj = 2, jpjm1 
    181             DO ji = fs_2, fs_jpim1   ! vector opt. 
    182                ! horizontal advective trends 
    183 #if defined key_zco 
    184                zbtr = e1e2tr(ji,jj) 
    185 #else 
    186                zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) 
    187 #endif 
    188                zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
    189                   &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
    190                ! add it to the general tracer trends 
    191                pta(ji,jj,jk) = pta(ji,jj,jk) + zta 
    192             END DO 
    193          END DO 
    194          !                                             ! =============== 
    195       END DO                                           !   End of slab 
    196       !                                                ! =============== 
    197  
    198       ! Horizontal trend used in tra_adv_ztvd subroutine 
    199       zltu(:,:,:) = pta(:,:,:) - ztrdt(:,:,:)  
    200  
    201       !  Save the horizontal advective trends for diagnostic 
    202       ! ----------------------------------------------------- 
    203       IF( l_trdtra ) THEN 
     131               zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) ) 
     132               zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) ) 
     133               ! UBS scheme 
     134               zwx(ji,jj,jk) = 0.5 * (  zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk)  ) 
     135               zwy(ji,jj,jk) = 0.5 * (  zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk)  ) 
     136            END DO 
     137         END DO 
     138      END DO 
     139 
     140      zltu(:,:,:) = pta(:,:,:)      ! store pta trends 
     141 
     142      ! Horizontal advective trends 
     143      DO jk = 1, jpkm1 
     144         DO jj = 2, jpjm1 
     145            DO ji = fs_2, fs_jpim1   ! vector opt. 
     146               zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     147               ! horizontal advective trends added to the general tracer trends 
     148               pta(ji,jj,jk) = pta(ji,jj,jk) - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
     149                  &                                    + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
     150            END DO 
     151         END DO 
     152      END DO 
     153 
     154      zltu(:,:,:) = pta(:,:,:) - zltu(:,:,:)    ! store the Horizontal advective trend (used in tra_adv_ztvd subroutine) 
     155 
     156  
     157      IF( l_trdtra ) THEN                       ! trend diagnostics (from advective fluxes) 
    204158         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn ) 
    205159         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn ) 
    206160      ENDIF 
    207  
    208       ! "Poleward" heat or salt transport  
     161      !                                         ! "Poleward" heat or salt transports 
    209162      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
    210          IF( lk_zco ) THEN 
    211             DO jk = 1, jpkm1 
    212                DO jj = 2, jpjm1 
    213                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    214                     zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 
    215                   END DO 
    216                END DO 
    217             END DO 
    218          ENDIF 
    219          IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
    220          IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) ) 
    221       ENDIF 
    222  
     163         IF( ktra == jp_tem )   pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
     164         IF( ktra == jp_sal )   pst_adv(:) = ptr_vj( zwy(:,:,:) ) 
     165      ENDIF 
     166      !                                         ! control print 
    223167      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - had: ', mask1=tmask, clinfo3=cdtype ) 
    224168 
     
    226170      ! II. Vertical advection 
    227171      ! ---------------------- 
    228       IF( l_trdtra )   ztrdt(:,:,:) = pta(:,:,:)          ! Save ta and sa trends 
     172      IF( l_trdtra )   zltv(:,:,:) = pta(:,:,:)          ! store pta if trend diag. 
    229173     
    230174      ! TVD scheme the vertical direction   
    231175      CALL tra_adv_ztvd( kt, pwn, zltu, ptb, ptn, pta ) 
    232176 
    233       IF( l_trdtra )   THEN         ! vertical advective trend diagnostics 
    234          DO jk = 1, jpkm1 
     177      IF( l_trdtra )   THEN                     ! vertical advective trend diagnostics 
     178         DO jk = 1, jpkm1                       ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 
    235179            DO jj = 2, jpjm1 
    236180               DO ji = fs_2, fs_jpim1   ! vector opt. 
    237                   z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) / fse3t(ji,jj,jk) 
    238                   ztrdt(ji,jj,jk) = pta(ji,jj,jk) - ztrdt(ji,jj,jk+1)  +  ptn(ji,jj,jk) * z_hdivn  
     181                  zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     182                  z_hdivn = (  pwn(ji,jj,jk) - pwn(ji,jj,jk+1)  ) * zbtr 
     183                  zltv(ji,jj,jk) = pta(ji,jj,jk) - zltv(ji,jj,jk+1)  +  ptn(ji,jj,jk) * z_hdivn  
    239184               END DO 
    240185            END DO 
    241186         END DO 
    242          CALL trd_tra( kt, ktra, jpt_trd_zad, cdtype, ptrd3d=ztrdt ) 
    243          ! 
    244       ENDIF 
    245       
     187         CALL trd_tra( kt, ktra, jpt_trd_zad, cdtype, ptrd3d=zltv ) 
     188      ENDIF 
     189      !                                         ! control print 
    246190      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - zad: ', mask1=tmask, clinfo3=cdtype ) 
    247191      ! 
     
    249193 
    250194 
    251    SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, ptb, ptn, pta ) 
     195   SUBROUTINE tra_adv_ztvd( kt, pwn, phtrd, ptb, ptn, pta ) 
    252196      !!---------------------------------------------------------------------- 
    253197      !!                  ***  ROUTINE tra_adv_ztvd  *** 
    254198      !!  
    255       !! **  Purpose :   Compute the now trend due to total advection of  
    256       !!       tracers and add it to the general trend of tracer equations 
     199      !! **  Purpose :   Compute the vertical advective trend of a tracer 
     200      !!               and add it to its general trend  
    257201      !! 
    258202      !! **  Method  :   TVD scheme, i.e. 2nd order centered scheme with 
    259203      !!       corrected flux (monotonic correction) 
    260       !!       note: - this advection scheme needs a leap-frog time scheme 
    261       !! 
    262       !! ** Action : - update (ta,sa) with the now advective tracer trends 
    263       !!             - save the trends in (ztrdt,ztrds) ('key_trdtra') 
    264       !!---------------------------------------------------------------------- 
    265       INTEGER , INTENT(in   )                         ::   kt      ! ocean time-step 
    266       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn     ! verical effective velocity 
    267       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   zttrd   ! lateral advective trends on T & S  
    268       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn        ! before and now tracer fields 
    269       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend  
     204      !! 
     205      !! ** Action : - update pta with the vertical advective trend 
     206      !!             - trend diagnostic (lk_trdtra=T) 
     207      !!---------------------------------------------------------------------- 
     208      INTEGER , INTENT(in   )                         ::   kt         ! ocean time-step 
     209      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwn        ! verical effective velocity 
     210      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   phtrd      ! lateral advective trends on T & S  
     211      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn   ! before and now tracer fields 
     212      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta        ! tracer trend  
    270213      !! 
    271214      INTEGER  ::   ji, jj, jk              ! dummy loop indices 
    272       REAL(wp) ::   z2dtt, zbtr, zew, z2    ! temporary scalar   
    273       REAL(wp) ::   ztak, zfp_wk, zfm_wk            !    "         " 
     215      REAL(wp) ::   z2dtt, zbtr, z2         ! temporary scalar   
     216      REAL(wp) ::   ztak, zfp_wk, zfm_wk    !    "         " 
    274217      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztw   ! temporary 3D workspace 
    275218      !!---------------------------------------------------------------------- 
     
    285228      ENDIF 
    286229 
    287       !  Bottom value : flux set to zero 
    288       ! -------------- 
    289       ztw(:,:,jpk) = 0.e0   ;   zti  (:,:,:) = 0.e0 
    290  
    291230 
    292231      !  upstream advection with initial mass fluxes & intermediate update 
    293       ! ------------------------------------------------------------------- 
    294       ! Surface value 
    295       IF( lk_dynspg_rl .OR. lk_vvl ) THEN                           ! rigid lid : flux set to zero 
    296          ztw(:,:,1) = 0.e0 
    297       ELSE                                                          ! free surface 
    298          ztw(:,:,1) = e1t(:,:) * e2t(:,:) * pwn(:,:,1) * ptb(:,:,1) 
    299       ENDIF 
    300  
    301       ! Interior value 
    302       DO jk = 2, jpkm1 
     232      ztw(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero 
     233      zti(:,:,jpk) = 0.e0 
     234      !                                                     ! Surface value :  
     235      IF( lk_dynspg_rl .OR. lk_vvl ) THEN   ;   ztw(:,:, 1 ) = 0.e0                      ! rigid lid or non-linear fs 
     236      ELSE                                  ;   ztw(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)   ! linear free surface  
     237      ENDIF 
     238      ! 
     239      DO jk = 2, jpkm1                                      ! interior values 
    303240         DO jj = 1, jpj 
    304241            DO ji = 1, jpi 
    305                zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 
    306                zfp_wk = zew + ABS( zew ) 
    307                zfm_wk = zew - ABS( zew ) 
    308                ztw(ji,jj,jk) = zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1) 
     242               zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     243               zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     244               ztw(ji,jj,jk) = 0.5 * (  zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1)  ) 
    309245            END DO 
    310246         END DO 
     
    318254               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    319255               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 
    320                pta(ji,jj,jk) =  pta(ji,jj,jk) + ztak 
    321                zti (ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * ( ztak + zttrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
    322             END DO 
    323          END DO 
    324       END DO 
    325  
    326       ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
    327       CALL lbc_lnk( zti, 'T', 1. ) 
     256               pta(ji,jj,jk) =   pta(ji,jj,jk) +           ztak 
     257               zti(ji,jj,jk) = ( ptb(ji,jj,jk) + z2dtt * ( ztak + phtrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
     258            END DO 
     259         END DO 
     260      END DO 
     261      ! 
     262      CALL lbc_lnk( zti, 'T', 1. )      ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
    328263 
    329264 
    330265      !  antidiffusive flux : high order minus low order 
    331       ! -------------------------------------------------       
    332266      ztw(:,:,1) = 0.e0       ! Surface value 
    333  
    334267      DO jk = 2, jpkm1        ! Interior value 
    335268         DO jj = 1, jpj 
    336269            DO ji = 1, jpi 
    337                zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 
    338                ztw(ji,jj,jk) = zew * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 
    339             END DO 
    340          END DO 
    341       END DO 
    342  
    343       !  monotonicity algorithm 
    344       ! ------------------------ 
    345       CALL nonosc_z( ptb, ztw, zti, z2 ) 
     270               ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 
     271            END DO 
     272         END DO 
     273      END DO 
     274      ! 
     275      CALL nonosc_z( ptb, ztw, zti, z2 )      !  monotonicity algorithm 
    346276 
    347277 
    348278      !  final trend with corrected fluxes 
    349       ! ----------------------------------- 
    350279      DO jk = 1, jpkm1 
    351280         DO jj = 2, jpjm1 
    352281            DO ji = fs_2, fs_jpim1   ! vector opt.   
    353282               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    354                ! k- vertical advective trends 
    355                ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 
    356                ! add them to the general tracer trends 
    357                pta(ji,jj,jk) = pta(ji,jj,jk) + ztak 
     283               ! k- vertical advective trends added to the general tracer trends 
     284               pta(ji,jj,jk) = pta(ji,jj,jk) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 
    358285            END DO 
    359286         END DO 
Note: See TracChangeset for help on using the changeset viewer.