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_ubs.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_ubs.F90

    r719 r786  
    22   !!============================================================================== 
    33   !!                       ***  MODULE  traadv_ubs  *** 
    4    !! Ocean active tracers:  horizontal & vertical advective trend 
     4   !! Ocean tracers:  horizontal & vertical advective trend 
    55   !!============================================================================== 
    6    !! History :  9.0  !  06-08  (L. Debreu, R. Benshila)  Original code 
     6   !! History :  1.0  !  06-08  (L. Debreu, R. Benshila)  Original code 
     7   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC 
    78   !!---------------------------------------------------------------------- 
    89 
     
    1112   !!                 advection trends using a third order biaised scheme   
    1213   !!---------------------------------------------------------------------- 
    13    USE oce             ! ocean dynamics and active tracers 
    1414   USE dom_oce         ! ocean space and time domain 
    1515   USE trdmod 
     
    3333#  include "vectopt_loop_substitute.h90" 
    3434   !!---------------------------------------------------------------------- 
    35    !!   OPA 9.0 , LOCEAN-IPSL (2006)  
    36    !! $Header$ 
     35   !! NEMO/OPA & TRP 2.4 , LOCEAN-IPSL (2008)  
     36   !! $Id:$  
    3737   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)  
    3838   !!---------------------------------------------------------------------- 
     
    4040CONTAINS 
    4141 
    42    SUBROUTINE tra_adv_ubs( kt, pun, pvn, pwn ) 
     42   SUBROUTINE tra_adv_ubs( kt, cdtype, ktra, pun, pvn, pwn,   & 
     43      &                                      ptb, ptn, pta ) 
    4344      !!---------------------------------------------------------------------- 
    4445      !!                  ***  ROUTINE tra_adv_ubs  *** 
     
    7071      !! 
    7172      !! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.  
    72       !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.  
    73       !!---------------------------------------------------------------------- 
    74       USE oce, ONLY :   zwx => ua   ! use ua as workspace 
    75       USE oce, ONLY :   zwy => va   ! use va as workspace 
    76       !! 
    77       INTEGER , INTENT(in)                         ::   kt             ! ocean time-step index 
    78       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pun   ! effective ocean velocity, u_component 
    79       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pvn   ! effective ocean velocity, v_component 
    80       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::  pwn   ! effective ocean velocity, w_component 
     73      !!             Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731–1741.  
     74      !!---------------------------------------------------------------------- 
     75      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index 
     76      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator) 
     77      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index 
     78      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! 3 ocean velocity components 
     79      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptb, ptn        ! before and now tracer fields 
     80      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend  
    8181      !! 
    8282      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    83       REAL(wp) ::   zta, zsa, zbtr, zcoef                  ! temporary scalars 
    84       REAL(wp) ::   zfui, zfp_ui, zfm_ui, zcenut, zcenus   !    "         " 
    85       REAL(wp) ::   zfvj, zfp_vj, zfm_vj, zcenvt, zcenvs   !    "         " 
    86       REAL(wp) ::   z_hdivn_x, z_hdivn_y, z_hdivn          !    "         " 
     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                           !    "         " 
    8787      REAL(wp), DIMENSION(jpi,jpj)     ::   zeeu, zeev     ! temporary 2D workspace 
    88       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz , zww                        ! temporary 3D workspace 
     88      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwx, zwy              ! temporary 3D workspace 
    8989      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztu , ztv , zltu , zltv, ztrdt   !    "              " 
    90       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsu , zsv , zlsu , zlsv, ztrds   !    "              " 
    9190      !!---------------------------------------------------------------------- 
    9291 
    9392      zltu(:,:,:) = 0.e0 
    9493      zltv(:,:,:) = 0.e0 
    95       zlsu(:,:,:) = 0.e0 
    96       zlsv(:,:,:) = 0.e0 
    9794 
    9895      IF( kt == nit000 ) THEN 
     
    104101      ENDIF 
    105102 
    106       ! Save ta and sa trends 
    107       ztrdt(:,:,:) = ta(:,:,:) 
    108       ztrds(:,:,:) = sa(:,:,:) 
     103      ! store pta trends 
     104      ztrdt(:,:,:) = pta(:,:,:) 
    109105 
    110106      zcoef = 1./6. 
     
    132128         DO jj = 1, jpjm1 
    133129            DO ji = 1, fs_jpim1   ! vector opt. 
    134                ztu(ji,jj,jk) = zeeu(ji,jj) * ( tb(ji+1,jj  ,jk) - tb(ji,jj,jk) ) 
    135                zsu(ji,jj,jk) = zeeu(ji,jj) * ( sb(ji+1,jj  ,jk) - sb(ji,jj,jk) ) 
    136                ztv(ji,jj,jk) = zeev(ji,jj) * ( tb(ji  ,jj+1,jk) - tb(ji,jj,jk) ) 
    137                zsv(ji,jj,jk) = zeev(ji,jj) * ( sb(ji  ,jj+1,jk) - sb(ji,jj,jk) ) 
     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) ) 
    138132            END DO 
    139133         END DO 
     
    145139#endif          
    146140               zltu(ji,jj,jk) = (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)  ) * zcoef 
    147                zlsu(ji,jj,jk) = (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk)  ) * zcoef 
    148141               zltv(ji,jj,jk) = (  ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) * zcoef 
    149                zlsv(ji,jj,jk) = (  zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  ) * zcoef 
    150142            END DO 
    151143         END DO 
     
    155147 
    156148      ! Lateral boundary conditions on the laplacian (zlt,zls)   (unchanged sgn) 
    157       CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zlsu, 'T', 1. ) 
    158       CALL lbc_lnk( zltv, 'T', 1. )   ;    CALL lbc_lnk( zlsv, 'T', 1. ) 
     149      CALL lbc_lnk( zltu, 'T', 1. )   ;    CALL lbc_lnk( zltv, 'T', 1. ) 
    159150 
    160151      !                                                ! =============== 
     
    178169               zfm_vj = zfvj - ABS( zfvj ) 
    179170               ! centered scheme 
    180                zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) ) 
    181                zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) ) 
    182                zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) ) 
    183                zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) ) 
     171               zcenut = zfui * ( ptn(ji,jj,jk) + ptn(ji+1,jj  ,jk) ) 
     172               zcenvt = zfvj * ( ptn(ji,jj,jk) + ptn(ji  ,jj+1,jk) ) 
    184173               ! mixed centered / upstream scheme 
    185174               zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) -zfm_ui * zltu(ji+1,jj,jk) 
    186175               zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) -zfm_vj * zltv(ji,jj+1,jk) 
    187                zww(ji,jj,jk) = zcenus - zfp_ui * zlsu(ji,jj,jk) -zfm_ui * zlsu(ji+1,jj,jk) 
    188                zwz(ji,jj,jk) = zcenvs - zfp_vj * zlsv(ji,jj,jk) -zfm_vj * zlsv(ji,jj+1,jk) 
    189176            END DO 
    190177         END DO 
     
    201188               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   & 
    202189                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  ) 
    203                zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   & 
    204                   &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  ) 
    205190               ! add it to the general tracer trends 
    206                ta(ji,jj,jk) = ta(ji,jj,jk) + zta 
    207                sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 
     191               pta(ji,jj,jk) = pta(ji,jj,jk) + zta 
    208192            END DO 
    209193         END DO 
     
    213197 
    214198      ! Horizontal trend used in tra_adv_ztvd subroutine 
    215       zltu(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)  
    216       zlsu(:,:,:) = sa(:,:,:) - ztrds(:,:,:)  
    217  
    218       ! 3. Save the horizontal advective trends for diagnostic 
    219       ! ------------------------------------------------------ 
    220       IF( l_trdtra )   THEN 
    221          ! Recompute the hoizontal advection zta & zsa trends computed  
    222          ! at the step 2. above in making the difference between the new  
    223          ! trends and the previous one ta()/sa - ztrdt()/ztrds() and add 
    224          ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends 
    225          ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0 
    226          ! 
    227          ! T/S ZONAL advection trends 
    228          DO jk = 1, jpkm1 
    229             DO jj = 2, jpjm1 
    230                DO ji = fs_2, fs_jpim1   ! vector opt. 
    231                   !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 
    232 #if defined key_zco 
    233                   zbtr = e1e2tr(ji,jj) 
    234                   z_hdivn_x = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)          & 
    235                      &         - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr 
    236 #else 
    237                   zbtr = e1e2tr(ji,jj) / fse3t(ji,jj,jk) 
    238                   z_hdivn_x = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          & 
    239                      &         - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr 
    240 #endif 
    241                   ztrdt(ji,jj,jk) = - ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_x 
    242                   ztrds(ji,jj,jk) = - ( zww(ji,jj,jk) - zww(ji-1,jj,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_x 
    243                END DO 
    244             END DO 
    245          END DO 
    246          CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)    ! save the trends 
    247          ! 
    248          ! T/S MERIDIONAL advection trends 
    249          DO jk = 1, jpkm1 
    250             DO jj = 2, jpjm1 
    251                DO ji = fs_2, fs_jpim1   ! vector opt. 
    252                   !-- Compute merid. divergence by splitting hdivn (see divcur.F90) 
    253 #if defined key_zco 
    254                   zbtr      = e1e2tr(ji,jj) 
    255                   z_hdivn_y = (  e1v(ji,  jj) * pvn(ji,jj  ,jk)          & 
    256                      &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr 
    257 #else 
    258                   zbtr      = e1e2tr(ji,jj) / fse3t(ji,jj,jk) 
    259                   z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          & 
    260                      &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 
    261 #endif 
    262                   ztrdt(ji,jj,jk) = - ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) ) * zbtr + tn(ji,jj,jk) * z_hdivn_y           
    263                   ztrds(ji,jj,jk) = - ( zwz(ji,jj,jk) - zwz(ji,jj-1,jk) ) * zbtr + sn(ji,jj,jk) * z_hdivn_y 
    264                END DO 
    265             END DO 
    266          END DO 
    267          CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)     ! save the trends 
    268          ! 
    269       ENDIF 
    270  
    271       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs had  - Ta: ', mask1=tmask,   & 
    272          &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    273  
    274       ! "zonal" mean advective heat and salt transport  
    275       IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
     199      zltu(:,:,:) = pta(:,:,:) - ztrdt(:,:,:)  
     200 
     201      !  Save the horizontal advective trends for diagnostic 
     202      ! ----------------------------------------------------- 
     203      IF( l_trdtra ) THEN 
     204         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptn ) 
     205         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptn ) 
     206      ENDIF 
     207 
     208      ! "Poleward" heat or salt transport  
     209      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
    276210         IF( lk_zco ) THEN 
    277211            DO jk = 1, jpkm1 
    278212               DO jj = 2, jpjm1 
    279213                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    280                      zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 
    281                      zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk) 
     214                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 
    282215                  END DO 
    283216               END DO 
    284217            END DO 
    285218         ENDIF 
    286          pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
    287          pst_adv(:) = ptr_vj( zwz(:,:,:) ) 
    288       ENDIF 
     219         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
     220         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) ) 
     221      ENDIF 
     222 
     223      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - had: ', mask1=tmask, clinfo3=cdtype ) 
     224 
    289225 
    290226      ! II. Vertical advection 
    291227      ! ---------------------- 
    292       IF( l_trdtra ) THEN          ! Save ta and sa trends 
    293          ztrdt(:,:,:) = ta(:,:,:) 
    294          ztrds(:,:,:) = sa(:,:,:) 
    295       ENDIF 
     228      IF( l_trdtra )   ztrdt(:,:,:) = pta(:,:,:)          ! Save ta and sa trends 
    296229     
    297230      ! TVD scheme the vertical direction   
    298       CALL tra_adv_ztvd(kt, pwn, zltu, zlsu) 
    299  
    300       IF( l_trdtra )   THEN         !  Save the final vertical advective trends 
     231      CALL tra_adv_ztvd( kt, pwn, zltu, ptb, ptn, pta ) 
     232 
     233      IF( l_trdtra )   THEN         !  vertical advective trend diagnostics 
    301234         DO jk = 1, jpkm1 
    302235            DO jj = 2, jpjm1 
    303236               DO ji = fs_2, fs_jpim1   ! vector opt. 
    304 #if defined key_zco 
    305                   zbtr      = e1e2tr(ji,jj) 
    306                   z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) 
    307                   z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) 
    308 #else 
    309                   zbtr      = e1e2tr(ji,jj) / fse3t(ji,jj,jk) 
    310                   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) 
    311                   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) 
    312 #endif 
    313                   z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr 
    314                   zbtr      = e1e2tr(ji,jj) / fse3t(ji,jj,jk) 
    315                   ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn 
    316                   ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn 
     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  
    317239               END DO 
    318240            END DO 
    319241         END DO 
    320          CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)   ! <<< ADD TO PREVIOUSLY COMPUTED 
     242         CALL trd_tra( kt, ktra, jpt_trd_zad, cdtype, ptrd3d=ztrdt ) 
    321243         ! 
    322244      ENDIF 
    323  
    324       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' ubs zad  - Ta: ', mask1=tmask,   & 
    325          &                       tab3d_2=sa, clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra') 
     245      
     246      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' ubs - zad: ', mask1=tmask, clinfo3=cdtype ) 
    326247      ! 
    327248   END SUBROUTINE tra_adv_ubs 
    328249 
    329250 
    330    SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, zstrd ) 
     251   SUBROUTINE tra_adv_ztvd( kt, pwn, zttrd, ptb, ptn, pta ) 
    331252      !!---------------------------------------------------------------------- 
    332253      !!                  ***  ROUTINE tra_adv_ztvd  *** 
     
    342263      !!             - save the trends in (ztrdt,ztrds) ('key_trdtra') 
    343264      !!---------------------------------------------------------------------- 
    344       INTEGER , INTENT(in)                         ::   kt             ! ocean time-step 
    345       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn            ! verical effective velocity 
    346       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   zttrd, zstrd   ! lateral advective trends on T & S  
     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  
    347270      !! 
    348271      INTEGER  ::   ji, jj, jk              ! dummy loop indices 
    349272      REAL(wp) ::   z2dtt, zbtr, zew, z2    ! temporary scalar   
    350       REAL(wp) ::   ztak, zfp_wk            !    "         " 
    351       REAL(wp) ::   zsak, zfm_wk            !    "         " 
     273      REAL(wp) ::   ztak, zfp_wk, zfm_wk            !    "         " 
    352274      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zti, ztw   ! temporary 3D workspace 
    353       REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zsi, zsw   !    "              " 
    354275      !!---------------------------------------------------------------------- 
    355276 
     
    366287      !  Bottom value : flux set to zero 
    367288      ! -------------- 
    368       ztw(:,:,jpk) = 0.e0   ;   zsw(:,:,jpk) = 0.e0 
    369       zti  (:,:,:) = 0.e0   ;   zsi  (:,:,:) = 0.e0 
     289      ztw(:,:,jpk) = 0.e0   ;   zti  (:,:,:) = 0.e0 
    370290 
    371291 
     
    375295      IF( lk_dynspg_rl .OR. lk_vvl ) THEN                           ! rigid lid : flux set to zero 
    376296         ztw(:,:,1) = 0.e0 
    377          zsw(:,:,1) = 0.e0 
    378       ELSE                                              ! free surface 
    379          DO jj = 1, jpj 
    380             DO ji = 1, jpi 
    381                zew = e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,1) 
    382                ztw(ji,jj,1) = zew * tb(ji,jj,1) 
    383                zsw(ji,jj,1) = zew * sb(ji,jj,1) 
    384             END DO 
    385          END DO 
     297      ELSE                                                          ! free surface 
     298         ztw(:,:,1) = e1t(:,:) * e2t(:,:) * pwn(:,:,1) * ptb(:,:,1) 
    386299      ENDIF 
    387300 
     
    393306               zfp_wk = zew + ABS( zew ) 
    394307               zfm_wk = zew - ABS( zew ) 
    395                ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1) 
    396                zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1) 
     308               ztw(ji,jj,jk) = zfp_wk * ptb(ji,jj,jk) + zfm_wk * ptb(ji,jj,jk-1) 
    397309            END DO 
    398310         END DO 
     
    406318               zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    407319               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 
    408                zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr 
    409                ta(ji,jj,jk) =  ta(ji,jj,jk) + ztak 
    410                sa(ji,jj,jk) =  sa(ji,jj,jk) + zsak  
    411                zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * ( ztak + zttrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
    412                zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * ( zsak + zstrd(ji,jj,jk) ) ) * tmask(ji,jj,jk) 
     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) 
    413322            END DO 
    414323         END DO 
     
    417326      ! Lateral boundary conditions on zti, zsi   (unchanged sign) 
    418327      CALL lbc_lnk( zti, 'T', 1. ) 
    419       CALL lbc_lnk( zsi, 'T', 1. ) 
    420328 
    421329 
    422330      !  antidiffusive flux : high order minus low order 
    423331      ! -------------------------------------------------       
    424       ! Surface value 
    425       ztw(:,:,1) = 0.e0   ;   zsw(:,:,1) = 0.e0 
    426  
    427       ! Interior value 
    428       DO jk = 2, jpkm1 
     332      ztw(:,:,1) = 0.e0       ! Surface value 
     333 
     334      DO jk = 2, jpkm1        ! Interior value 
    429335         DO jj = 1, jpj 
    430336            DO ji = 1, jpi 
    431337               zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * pwn(ji,jj,jk) 
    432                ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 
    433                zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk) 
     338               ztw(ji,jj,jk) = zew * ( ptn(ji,jj,jk) + ptn(ji,jj,jk-1) ) - ztw(ji,jj,jk) 
    434339            END DO 
    435340         END DO 
     
    438343      !  monotonicity algorithm 
    439344      ! ------------------------ 
    440       CALL nonosc_z( tb, ztw, zti, z2 ) 
    441       CALL nonosc_z( sb, zsw, zsi, z2 ) 
     345      CALL nonosc_z( ptb, ztw, zti, z2 ) 
    442346 
    443347 
     
    450354               ! k- vertical advective trends 
    451355               ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 
    452                zsak = - ( zsw(ji,jj,jk) - zsw(ji,jj,jk+1) ) * zbtr 
    453356               ! add them to the general tracer trends 
    454                ta(ji,jj,jk) = ta(ji,jj,jk) + ztak 
    455                sa(ji,jj,jk) = sa(ji,jj,jk) + zsak 
     357               pta(ji,jj,jk) = pta(ji,jj,jk) + ztak 
    456358            END DO 
    457359         END DO 
Note: See TracChangeset for help on using the changeset viewer.