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 786 for branches/dev_001_GM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 – NEMO

Ignore:
Timestamp:
2008-01-10T18:11:23+01:00 (16 years ago)
Author:
gm
Message:

dev_001_GM - merge TRC-TRA on OPA only, trabbl & zpshde not done and trdmld not OK - compilation OK

File:
1 edited

Legend:

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

    r719 r786  
    66   !! History :  8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv  
    77   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
    8    !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization 
    9    !!            " "  !  06-04  (R. Benshila, G. Madec) Step reorganization 
    10    !!---------------------------------------------------------------------- 
    11  
    12    !!---------------------------------------------------------------------- 
    13    !!   tra_adv_cen2 : update the tracer trend with the horizontal and 
    14    !!                  vertical advection trends using a seconder order 
    15    !!                  centered scheme. (k-j-i loops) 
    16    !!---------------------------------------------------------------------- 
    17    USE oce             ! ocean dynamics and active tracers 
     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 
     12   !!---------------------------------------------------------------------- 
     13 
     14   !!---------------------------------------------------------------------- 
     15   !!   tra_adv_cen2 : update the tracer trend with the horizontal and vertical 
     16   !!                  advection trends using a 2nd order centered scheme 
     17   !!---------------------------------------------------------------------- 
     18   USE oce, ONLY: tn   ! now ocean temperature 
    1819   USE dom_oce         ! ocean space and time domain 
    1920   USE trdmod          ! ocean active tracers trends  
    2021   USE trdmod_oce      ! ocean variables trends 
    2122   USE flxrnf          ! 
    22    USE trabbl          ! advective term in the BBL 
    2323   USE ocfzpt          ! 
    2424   USE lib_mpp 
     
    4040#  include "vectopt_loop_substitute.h90" 
    4141   !!---------------------------------------------------------------------- 
    42    !!   OPA 9.0 , LOCEAN-IPSL (2005)  
    43    !! $Header$  
     42   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)  
     43   !! $Id:$  
    4444   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4545   !!---------------------------------------------------------------------- 
     
    4747CONTAINS 
    4848 
    49    SUBROUTINE tra_adv_cen2( kt, pun, pvn, pwn ) 
     49   SUBROUTINE tra_adv_cen2( kt, cdtype, ktra, pun, pvn, pwn,   & 
     50      &                                       ptb, ptn, pta ) 
    5051      !!---------------------------------------------------------------------- 
    5152      !!                  ***  ROUTINE tra_adv_cen2  *** 
    5253      !!                  
    53       !! ** Purpose :   Compute the now trend due to the advection of tracers 
    54       !!      and add it to the general trend of passive tracer equations. 
     54      !! ** Purpose :   Compute the now trend due to the advection of a tracer 
     55      !!      and add it to the corresponding general trend of tracer equations. 
    5556      !! 
    5657      !! ** Method  :   The advection is evaluated by a second order centered 
     
    6364      !!         Part I : horizontal advection 
    6465      !!       * centered flux: 
    65       !!               zcenu = e2u*e3u  un  mi(tn) 
    66       !!               zcenv = e1v*e3v  vn  mj(tn) 
     66      !!               zcenu = e2u*e3u  un  mi(ptn) 
     67      !!               zcenv = e1v*e3v  vn  mj(ptn) 
    6768      !!       * upstream flux: 
    68       !!               zupsu = e2u*e3u  un  (tb(i) or tb(i-1) ) [un>0 or <0] 
     69      !!               zupsu = e2u*e3u  un  (ptb(i) or ptb(i-1) ) [un>0 or <0] 
    6970      !!               zupsv = e1v*e3v  vn  (tb(j) or tb(j-1) ) [vn>0 or <0] 
    7071      !!       * mixed upstream / centered horizontal advection scheme 
     
    7576      !!       * horizontal advective trend (divergence of the fluxes) 
    7677      !!               zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] } 
    77       !!       * Add this trend now to the general trend of tracer (ta,sa): 
    78       !!              (ta,sa) = (ta,sa) + ( zta , zsa ) 
    79       !!       * trend diagnostic ('key_trdtra' defined): the trend is 
     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 
    8081      !!      saved for diagnostics. The trends saved is expressed as 
    8182      !!      Uh.gradh(T), i.e. 
    82       !!                     save trend = zta + tn divn 
     83      !!                     save trend = zta + ptn divn 
    8384      !!         In addition, the advective trend in the two horizontal direc- 
    84       !!      tion is also re-computed as Uh gradh(T). Indeed hadt+tn divn is 
     85      !!      tion is also re-computed as Uh gradh(T). Indeed hadt+ptn divn is 
    8586      !!      equal to (in s-coordinates, and similarly in z-coord.): 
    86       !!         zta+tn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[tn] ) 
    87       !!                                      +mj-1( e1v*e3v  vn  mj[tn] )  } 
     87      !!         zta+ptn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[ptn] ) 
     88      !!                                       +mj-1( e1v*e3v  vn  mj[ptn] )  } 
    8889      !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so 
    8990      !!      they vanish from the expression of the flux and divergence. 
    9091      !! 
    9192      !!         Part II : vertical advection 
    92       !!      For temperature (idem for salinity) the advective trend is com- 
    93       !!      puted as follows : 
    94       !!            zta = 1/e3t dk+1[ zwz ] 
    95       !!      where the vertical advective flux, zwz, is given by : 
    96       !!            zwz = zcofk * zupst + (1-zcofk) * zcent 
     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 
    9797      !!      with 
    9898      !!        zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0] 
    99       !!        zcenu = centered flux = wn * mk(tn) 
     99      !!        zcenu = centered flux = wn * mk(ptn) 
    100100      !!         The surface boundary condition is : 
    101101      !!      rigid-lid (lk_dynspg_frd = T) : zero advective flux 
    102       !!      free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * tn(:,:,1) 
     102      !!      free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * ptn(:,:,1) 
    103103      !!         Add this trend now to the general trend of tracer (ta,sa): 
    104104      !!            (ta,sa) = (ta,sa) + ( zta , zsa ) 
    105105      !!         Trend diagnostic ('key_trdtra' defined): the trend is 
    106106      !!      saved for diagnostics. The trends saved is expressed as : 
    107       !!             save trend =  w.gradz(T) = zta - tn divn. 
     107      !!             save trend =  w.gradz(T) = zta - ptn divn. 
    108108      !! 
    109109      !! ** Action :  - update (ta,sa) with the now advective tracer trends 
    110       !!              - save trends in (ztrdt,ztrds) ('key_trdtra') 
     110      !!              - trend diagnostics (lk_trdtra=T) 
    111111      !!---------------------------------------------------------------------- 
    112       USE oce, ONLY :   zwx => ua   ! use ua as workspace 
    113       USE oce, ONLY :   zwy => va   ! use va as workspace 
    114       !! 
    115       INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index 
    116       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component 
    117       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component 
    118       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component 
    119       !! 
    120       INTEGER  ::   ji, jj, jk                           ! dummy loop indices 
    121       REAL(wp) ::                           & 
    122          zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars 
    123          zhw, ze3tr, zcofi, zcofj,          &  !    "         " 
    124          zupsut, zupsvt, zupsus, zupsvs,    &  !    "         " 
    125          zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         " 
    126          zcofk, zupst, zupss, zcent,        &  !    "         " 
    127          zcens, zfp_w, zfm_w,               &  !    "         " 
    128          zcenut, zcenvt, zcenus, zcenvs,    &  !    "         " 
    129          z_hdivn_x, z_hdivn_y, z_hdivn 
    130       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind   ! 3D workspace  
    131       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww, ztrds         !  "      " 
     112      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index 
     113      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator) 
     114      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index 
     115      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! 3 ocean velocity components 
     116      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn        ! before and now tracer fields 
     117      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend  
     118      !! 
     119      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  
    132125      !!---------------------------------------------------------------------- 
    133126 
     
    135128         IF(lwp) WRITE(numout,*) 
    136129         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme' 
    137          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case' 
    138          IF(lwp) WRITE(numout,*) 
     130         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    139131         !  
    140132         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) ) 
     
    158150      END DO 
    159151 
    160  
    161       !  Horizontal advective fluxes 
    162       ! ----------------------------- 
     152      ! I. Horizontal advection 
     153      ! ----------------------- 
     154 
    163155      !                                                ! =============== 
    164156      DO jk = 1, jpkm1                                 ! Horizontal slab 
    165157         !                                             ! =============== 
    166          DO jj = 1, jpjm1 
     158         ! 
     159         DO jj = 1, jpjm1                              !  Horizontal advective fluxes 
    167160            DO ji = 1, fs_jpim1   ! vector opt. 
    168161               ! upstream indicator 
     
    182175               zfm_ui = zfui - ABS( zfui ) 
    183176               zfm_vj = zfvj - ABS( zfvj ) 
    184                zupsut = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj  ,jk) 
    185                zupsvt = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji  ,jj+1,jk) 
    186                zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk) 
    187                zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk) 
     177               zupsut = zfp_ui * ptb(ji,jj,jk) + zfm_ui * ptb(ji+1,jj  ,jk) 
     178               zupsvt = zfp_vj * ptb(ji,jj,jk) + zfm_vj * ptb(ji  ,jj+1,jk) 
    188179               ! centered scheme 
    189                zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) ) 
    190                zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) ) 
    191                zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) ) 
    192                zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) ) 
     180               zcenut = zfui * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) ) 
     181               zcenvt = zfvj * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) ) 
    193182               ! mixed centered / upstream scheme 
    194183               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut 
    195184               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt 
    196                zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus 
    197                zwz(ji,jj,jk) = zcofj * zupsvs + (1.-zcofj) * zcenvs 
    198             END DO 
    199          END DO 
    200  
    201          !  Tracer flux divergence at t-point added to the general trend 
    202          ! -------------------------------------------------------------- 
    203          DO jj = 2, jpjm1 
     185            END DO 
     186         END DO 
     187 
     188         DO jj = 2, jpjm1                              !  horizontal tracer flux divergence added to the general trend 
    204189            DO ji = fs_2, fs_jpim1   ! vector opt. 
    205190#if defined key_zco 
     
    211196               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
    212197                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
    213                zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   & 
    214                   &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  ) 
    215198               ! add it to the general tracer trends 
    216                ta(ji,jj,jk) = ta(ji,jj,jk) + zta 
    217                sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 
     199               pta(ji,jj,jk) = pta(ji,jj,jk) + zta 
    218200            END DO 
    219201         END DO 
     
    225207      ! ----------------------------------------------------- 
    226208      IF( l_trdtra ) THEN 
    227          ! T/S ZONAL advection trends 
    228          ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0 
    229          ! 
    230          DO jk = 1, jpkm1 
    231             DO jj = 2, jpjm1 
    232                DO ji = fs_2, fs_jpim1   ! vector opt. 
    233                   !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 
    234                   !   N.B. This computation is not valid along OBCs (if any) 
    235 #if defined key_zco 
    236                   zbtr      = btr2(ji,jj)  
    237                   z_hdivn_x = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)                              & 
    238                      &         - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr 
    239 #else 
    240                   zbtr      = btr2(ji,jj) / fse3t(ji,jj,jk) 
    241                   z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          & 
    242                      &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr 
    243 #endif 
    244                   ztrdt(ji,jj,jk) = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x 
    245                   ztrds(ji,jj,jk) = - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj,jk) ) + sn(ji,jj,jk) * z_hdivn_x 
    246                END DO 
    247             END DO 
    248          END DO 
    249          CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt) 
    250          ! 
    251          ! T/S MERIDIONAL advection trends 
    252          DO jk = 1, jpkm1 
    253             DO jj = 2, jpjm1 
    254                DO ji = fs_2, fs_jpim1   ! vector opt. 
    255                   !-- Compute merid. divergence by splitting hdivn (see divcur.F90) 
    256                   !   N.B. This computation is not valid along OBCs (if any) 
    257 #if defined key_zco 
    258                   zbtr      = btr2(ji,jj)  
    259                   z_hdivn_y = (  e1v(ji,jj  ) * pvn(ji,jj  ,jk)                              & 
    260                      &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr 
    261 #else 
    262                   zbtr      = btr2(ji,jj) / fse3t(ji,jj,jk) 
    263                   z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          & 
    264                      &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 
    265 #endif 
    266                   ztrdt(ji,jj,jk) = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y           
    267                   ztrds(ji,jj,jk) = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y 
    268                END DO 
    269             END DO 
    270          END DO 
    271          CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) 
    272          ! 
    273          ! Save the horizontal up-to-date ta/sa trends 
    274          ztrdt(:,:,:) = ta(:,:,:)  
    275          ztrds(:,:,:) = sa(:,:,:) 
    276       ENDIF 
    277  
    278       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 had  - Ta: ', mask1=tmask, & 
    279          &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    280  
    281       ! 4. "zonal" mean advective heat and salt transport  
    282       ! ------------------------------------------------- 
    283  
    284       IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) 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 
    285219         IF( lk_zco ) THEN 
    286220            DO jk = 1, jpkm1 
     
    288222                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    289223                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 
    290                     zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk) 
    291224                  END DO 
    292225               END DO 
    293226            END DO 
    294227         ENDIF 
    295          pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
    296          pst_adv(:) = ptr_vj( zwz(:,:,:) ) 
    297       ENDIF 
     228         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
     229         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) ) 
     230      ENDIF 
     231 
    298232 
    299233      ! II. Vertical advection 
    300234      ! ---------------------- 
    301235 
    302       ! Bottom value : flux set to zero 
    303       zwx(:,:,jpk) = 0.e0     ;    zwy(:,:,jpk) = 0.e0 
    304  
    305       ! Surface value 
    306       IF( lk_dynspg_rl .OR. lk_vvl ) THEN 
    307          ! rigid lid or variable volume: flux set to zero 
    308          zwx(:,:, 1 ) = 0.e0    ;    zwy(:,:, 1 ) = 0.e0 
    309       ELSE 
    310          ! free surface 
    311          zwx(:,:, 1 ) = pwn(:,:,1) * tn(:,:,1) 
    312          zwy(:,:, 1 ) = pwn(:,:,1) * sn(:,:,1) 
    313       ENDIF 
    314  
    315       ! 1. Vertical advective fluxes 
    316       ! ---------------------------- 
    317       ! Second order centered tracer flux at w-point 
    318       DO jk = 2, jpk 
     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) 
    319245         DO jj = 2, jpjm1 
    320246            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    326252               zfp_w = zhw + ABS( zhw ) 
    327253               zfm_w = zhw - ABS( zhw ) 
    328                zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1) 
    329                zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1) 
     254               zupst = zfp_w * ptb(ji,jj,jk) + zfm_w * ptb(ji,jj,jk-1) 
    330255               ! centered scheme 
    331                zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) 
    332                zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) 
     256               zcent = zhw * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) 
    333257               ! mixed centered / upstream scheme 
    334258               zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent 
    335                zwy(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens 
    336259            END DO 
    337260         END DO 
    338261      END DO 
    339262 
    340       ! 2. Tracer flux divergence at t-point added to the general trend 
    341       ! ------------------------- 
    342       DO jk = 1, jpkm1 
     263      DO jk = 1, jpkm1                                 ! Tracer flux divergence at t-point added to the general trend 
    343264         DO jj = 2, jpjm1 
    344265            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    346267               ! vertical advective trends 
    347268               zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
    348                zsa = - ze3tr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) 
    349269               ! add it to the general tracer trends 
    350                ta(ji,jj,jk) =  ta(ji,jj,jk) + zta 
    351                sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa 
     270               pta(ji,jj,jk) =  pta(ji,jj,jk) + zta 
    352271            END DO 
    353272         END DO 
     
    356275      ! 3. Save the vertical advective trends for diagnostic 
    357276      ! ---------------------------------------------------- 
    358       IF( l_trdtra )   THEN 
    359          ! Recompute the vertical advection zta & zsa trends computed  
    360          ! at the step 2. above in making the difference between the new  
    361          ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract 
    362          ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends 
    363  
    364          DO jk = 1, jpkm1 
    365             DO jj = 2, jpjm1 
    366                DO ji = fs_2, fs_jpim1   ! vector opt. 
    367 #if defined key_zco 
    368                   zbtr      = btr2(ji,jj)  
    369                   z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) 
    370                   z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) 
    371 #else 
    372                   zbtr      = btr2(ji,jj) / fse3t(ji,jj,jk) 
    373                   z_hdivn_x = e2u(ji,jj)*fse3u(ji,jj,jk)*pun(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk)*pun(ji-1,jj,jk) 
    374                   z_hdivn_y = e1v(ji,jj)*fse3v(ji,jj,jk)*pvn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk)*pvn(ji,jj-1,jk) 
    375 #endif 
    376                   z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr 
    377                   ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn  
    378                   ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn 
    379                END DO 
    380             END DO 
    381          END DO 
    382          CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) 
    383       ENDIF 
    384  
    385       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 zad  - Ta: ', mask1=tmask, & 
    386          &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     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 ) 
    387280      ! 
    388281   END SUBROUTINE tra_adv_cen2 
Note: See TracChangeset for help on using the changeset viewer.