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

    r719 r786  
    66   !! History :       !  06-00  (A.Estublier)  for passive tracers 
    77   !!                 !  01-08  (E.Durand, G.Madec)  adapted for T & S 
    8    !!            8.5  !  02-06  (G. Madec)  F90: Free form and module 
     8   !!   NEMO     1.0  !  02-06  (G. Madec)  F90: Free form and module 
     9   !!            2.4  !  08-01  (G. Madec) Merge TRA-TRC 
    910   !!---------------------------------------------------------------------- 
    1011 
     
    1314   !!                   and vertical advection trends using MUSCL scheme 
    1415   !!---------------------------------------------------------------------- 
    15    USE oce             ! ocean dynamics and active tracers 
    1616   USE dom_oce         ! ocean space and time domain 
    1717   USE trdmod          ! ocean active tracers trends  
     
    3434#  include "vectopt_loop_substitute.h90" 
    3535   !!---------------------------------------------------------------------- 
    36    !!   OPA 9.0 , LOCEAN-IPSL (2006)  
    37    !! $Header$  
     36   !! NEMO/OPA 2.4 , LOCEAN-IPSL (2008)  
     37   !! $Id:$  
    3838   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    3939   !!---------------------------------------------------------------------- 
     
    4141CONTAINS 
    4242 
    43    SUBROUTINE tra_adv_muscl( kt, pun, pvn, pwn ) 
     43   SUBROUTINE tra_adv_muscl( kt, cdtype, ktra, pun, pvn, pwn,   & 
     44      &                                        ptb     , pta ) 
    4445      !!---------------------------------------------------------------------- 
    4546      !!                    ***  ROUTINE tra_adv_muscl  *** 
     
    5152      !! ** Method  : MUSCL scheme plus centered scheme at ocean boundaries 
    5253      !! 
    53       !! ** Action  : - update (ta,sa) with the now advective tracer trends 
     54      !! ** Action  : - update (pta,sa) with the now advective tracer trends 
    5455      !!              - save trends in (ztrdt,ztrds) ('key_trdtra') 
    5556      !! 
     
    5758      !!              IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa) 
    5859      !!---------------------------------------------------------------------- 
    59       USE oce, ONLY :   ztrdt => ua   ! use ua as workspace 
    60       USE oce, ONLY :   ztrds => va   ! use va as workspace 
    61       !! 
    62       INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index 
    63       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component 
    64       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component 
    65       REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component 
    66       !! 
    67       INTEGER ::   ji, jj, jk   ! dummy loop indices 
    68       REAL(wp) ::   & 
    69          zu, zv, zw, zeu, zev,           &   
    70          zew, zbtr, zstep,               & 
    71          z0u, z0v, z0w,                  & 
    72          zzt1, zzt2, zalpha,             & 
    73          zzs1, zzs2, z2,                 & 
    74          zta, zsa,                       & 
    75         z_hdivn_x, z_hdivn_y, z_hdivn 
    76       REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zt1, zt2, ztp1, ztp2   ! 3D workspace 
    77       REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zs1, zs2, zsp1, zsp2   !  "      " 
     60      INTEGER         , INTENT(in   )                         ::   kt              ! ocean time-step index 
     61      CHARACTER(len=3), INTENT(in   )                         ::   cdtype          ! =TRA or TRC (tracer indicator) 
     62      INTEGER         , INTENT(in   )                         ::   ktra            ! tracer index 
     63      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! 3 ocean velocity components 
     64      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptb             ! before tracer fields 
     65      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pta             ! tracer trend  
     66      !! 
     67      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 
     73      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zwx, zwy, zslpx, zslpy   ! 3D workspace 
    7874      !!---------------------------------------------------------------------- 
    7975 
     
    9591         DO jj = 1, jpjm1       
    9692            DO ji = 1, fs_jpim1   ! vector opt. 
    97                zt1(ji,jj,jk) = umask(ji,jj,jk) * ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) 
    98                zs1(ji,jj,jk) = umask(ji,jj,jk) * ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) 
    99                zt2(ji,jj,jk) = vmask(ji,jj,jk) * ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) 
    100                zs2(ji,jj,jk) = vmask(ji,jj,jk) * ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) 
     93               zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk) - ptb(ji,jj,jk) ) 
     94               zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk) - ptb(ji,jj,jk) ) 
    10195            END DO 
    10296         END DO 
    10397      END DO 
    10498      ! bottom values 
    105       zt1(:,:,jpk) = 0.e0    ;    zt2(:,:,jpk) = 0.e0 
    106       zs1(:,:,jpk) = 0.e0    ;    zs2(:,:,jpk) = 0.e0 
    107  
    108       ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign) 
    109       CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. ) 
    110       CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. ) 
     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. ) 
    111103 
    112104      ! Slopes 
     
    115107         DO jj = 2, jpj 
    116108            DO ji = fs_2, jpi   ! vector opt. 
    117                ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji-1,jj  ,jk) )   & 
    118                   &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji-1,jj  ,jk) ) ) 
    119                zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji-1,jj  ,jk) )   & 
    120                   &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji-1,jj  ,jk) ) ) 
    121                ztp2(ji,jj,jk) =                    ( zt2(ji,jj,jk) + zt2(ji  ,jj-1,jk) )   & 
    122                   &           * ( 0.25 + SIGN( 0.25, zt2(ji,jj,jk) * zt2(ji  ,jj-1,jk) ) ) 
    123                zsp2(ji,jj,jk) =                    ( zs2(ji,jj,jk) + zs2(ji  ,jj-1,jk) )   & 
    124                   &           * ( 0.25 + SIGN( 0.25, zs2(ji,jj,jk) * zs2(ji  ,jj-1,jk) ) ) 
     109               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) ) ) 
     111               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) ) ) 
    125113            END DO 
    126114         END DO 
    127115      END DO 
    128116      ! bottom values 
    129       ztp1(:,:,jpk) = 0.e0    ;    ztp2(:,:,jpk) = 0.e0 
    130       zsp1(:,:,jpk) = 0.e0    ;    zsp2(:,:,jpk) = 0.e0 
     117      zslpx(:,:,jpk) = 0.e0    ;    zslpy(:,:,jpk) = 0.e0 
    131118 
    132119      ! Slopes limitation 
     
    134121         DO jj = 2, jpj 
    135122            DO ji = fs_2, jpi   ! vector opt. 
    136                ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   & 
    137                   &           * MIN(    ABS( ztp1(ji  ,jj,jk) ),   & 
    138                   &                  2.*ABS( zt1 (ji-1,jj,jk) ),   & 
    139                   &                  2.*ABS( zt1 (ji  ,jj,jk) ) ) 
    140                zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   & 
    141                   &           * MIN(    ABS( zsp1(ji  ,jj,jk) ),   & 
    142                   &                  2.*ABS( zs1 (ji-1,jj,jk) ),   & 
    143                   &                  2.*ABS( zs1 (ji  ,jj,jk) ) ) 
    144                ztp2(ji,jj,jk) = SIGN( 1., ztp2(ji,jj,jk) )   & 
    145                   &           * MIN(    ABS( ztp2(ji,jj  ,jk) ),   & 
    146                   &                  2.*ABS( zt2 (ji,jj-1,jk) ),   & 
    147                   &                  2.*ABS( zt2 (ji,jj  ,jk) ) ) 
    148                zsp2(ji,jj,jk) = SIGN( 1., zsp2(ji,jj,jk) )   & 
    149                   &           * MIN(    ABS( zsp2(ji,jj  ,jk) ),   & 
    150                   &                  2.*ABS( zs2 (ji,jj-1,jk) ),   & 
    151                   &                  2.*ABS( zs2 (ji,jj  ,jk) ) ) 
     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) ) ) 
    152131            END DO 
    153132         END DO 
     
    172151               zalpha = 0.5 - z0u 
    173152               zu  = z0u - 0.5 * pun(ji,jj,jk) * zstep / e1u(ji,jj) 
    174                zzt1 = tb(ji+1,jj,jk) + zu*ztp1(ji+1,jj,jk) 
    175                zzt2 = tb(ji  ,jj,jk) + zu*ztp1(ji  ,jj,jk) 
    176                zzs1 = sb(ji+1,jj,jk) + zu*zsp1(ji+1,jj,jk) 
    177                zzs2 = sb(ji  ,jj,jk) + zu*zsp1(ji  ,jj,jk) 
    178                zt1(ji,jj,jk) = zeu * ( zalpha * zzt1 + (1.-zalpha) * zzt2 ) 
    179                zs1(ji,jj,jk) = zeu * ( zalpha * zzs1 + (1.-zalpha) * zzs2 ) 
     153               zzwx = ptb(ji+1,jj,jk) + zu * zslpx(ji+1,jj,jk) 
     154               zzwy = ptb(ji  ,jj,jk) + zu * zslpx(ji  ,jj,jk) 
     155               zwx(ji,jj,jk) = zeu * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    180156               ! 
    181157               z0v = SIGN( 0.5, pvn(ji,jj,jk) )             
    182158               zalpha = 0.5 - z0v 
    183159               zv  = z0v - 0.5 * pvn(ji,jj,jk) * zstep / e2v(ji,jj) 
    184                zzt1 = tb(ji,jj+1,jk) + zv*ztp2(ji,jj+1,jk) 
    185                zzt2 = tb(ji,jj  ,jk) + zv*ztp2(ji,jj  ,jk) 
    186                zzs1 = sb(ji,jj+1,jk) + zv*zsp2(ji,jj+1,jk) 
    187                zzs2 = sb(ji,jj  ,jk) + zv*zsp2(ji,jj  ,jk) 
    188                zt2(ji,jj,jk) = zev * ( zalpha * zzt1 + (1.-zalpha) * zzt2 ) 
    189                zs2(ji,jj,jk) = zev * ( zalpha * zzs1 + (1.-zalpha) * zzs2 ) 
    190             END DO 
    191          END DO 
    192       END DO 
    193  
    194       ! lateral boundary conditions on zt1, zt2 ; zs1, zs2   (changed sign) 
    195       CALL lbc_lnk( zt1, 'U', -1. )   ;   CALL lbc_lnk( zs1, 'U', -1. )  
    196       CALL lbc_lnk( zt2, 'V', -1. )   ;   CALL lbc_lnk( zs2, 'V', -1. ) 
     160               zzwx = ptb(ji,jj+1,jk) + zv * zslpy(ji,jj+1,jk) 
     161               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) 
     169      CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. ) 
    197170 
    198171      ! Tracer flux divergence at t-point added to the general trend 
     
    201174            DO ji = fs_2, fs_jpim1   ! vector opt. 
    202175#if defined key_zco 
    203                zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj) ) 
     176               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 
    204177#else 
    205                zbtr = 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) ) 
     178               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    206179#endif 
    207180               ! horizontal advective trends 
    208                zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj  ,jk  )   & 
    209                   &           + zt2(ji,jj,jk) - zt2(ji  ,jj-1,jk  ) ) 
    210                zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj  ,jk  )   & 
    211                   &           + zs2(ji,jj,jk) - zs2(ji  ,jj-1,jk  ) )  
     181               zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     182                  &           + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  ) )  
    212183               ! add it to the general tracer trends 
    213                ta(ji,jj,jk) = ta(ji,jj,jk) + zta 
    214                sa(ji,jj,jk) = sa(ji,jj,jk) + zsa 
     184               pta(ji,jj,jk) = pta(ji,jj,jk) + zta 
    215185            END DO 
    216186        END DO 
    217187      END DO         
    218188 
    219       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl had  - Ta: ', mask1=tmask ,  & 
    220          &                       tab3d_2=sa, clinfo2=             ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     189      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - had: ', mask1=tmask, clinfo3=cdtype ) 
    221190 
    222191      ! Save the horizontal advective trends for diagnostics 
    223192      IF( l_trdtra ) THEN 
    224          ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0 
    225          ! 
    226          ! T/S ZONAL advection trends 
    227          DO jk = 1, jpkm1 
    228             DO jj = 2, jpjm1 
    229                DO ji = fs_2, fs_jpim1   ! vector opt. 
    230                   !-- Compute zonal divergence by splitting hdivn (see divcur.F90) 
    231                   !   N.B. This computation is not valid along OBCs (if any) 
    232 #if defined key_zco 
    233                   zbtr      = 1. / ( e1t(ji,jj) * e2t(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      = 1. / ( e1t(ji,jj) * e2t(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) = - zbtr * ( zt1(ji,jj,jk) - zt1(ji-1,jj,jk) ) + tn(ji,jj,jk) * z_hdivn_x 
    242                   ztrds(ji,jj,jk) = - zbtr * ( zs1(ji,jj,jk) - zs1(ji-1,jj,jk) ) + 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) 
    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                   !   N.B. This computation is not valid along OBCs (if any) 
    254 #if defined key_zco 
    255                   zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 
    256                   z_hdivn_y = (  e1v(ji,jj  ) * pvn(ji,jj  ,jk)                              & 
    257                      &         - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr 
    258 #else 
    259                   zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    260                   z_hdivn_y = (  e1v(ji,  jj) * fse3v(ji,jj  ,jk) * pvn(ji,jj  ,jk)          & 
    261                      &         - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr 
    262 #endif 
    263                   ztrdt(ji,jj,jk) = - zbtr * ( zt2(ji,jj,jk) - zt2(ji,jj-1,jk) ) + tn(ji,jj,jk) * z_hdivn_y           
    264                   ztrds(ji,jj,jk) = - zbtr * ( zs2(ji,jj,jk) - zs2(ji,jj-1,jk) ) + sn(ji,jj,jk) * z_hdivn_y 
    265                END DO 
    266             END DO 
    267          END DO 
    268          CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt) 
    269  
    270          ! Save the up-to-date ta and sa trends 
    271          ztrdt(:,:,:) = ta(:,:,:)  
    272          ztrds(:,:,:) = sa(:,:,:)  
    273          ! 
    274       ENDIF 
    275  
    276       ! "zonal" mean advective heat and salt transport 
    277       IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
     193         CALL trd_tra_adv( kt, ktra, jpt_trd_xad, cdtype, zwx, pun, ptb ) 
     194         CALL trd_tra_adv( kt, ktra, jpt_trd_yad, cdtype, zwy, pvn, ptb ) 
     195      ENDIF 
     196 
     197      IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN 
    278198         IF( lk_zco ) THEN 
    279199            DO jk = 1, jpkm1 
    280200               DO jj = 2, jpjm1 
    281201                  DO ji = fs_2, fs_jpim1   ! vector opt. 
    282                     zt2(ji,jj,jk) = zt2(ji,jj,jk) * fse3v(ji,jj,jk) 
    283                     zs2(ji,jj,jk) = zs2(ji,jj,jk) * fse3v(ji,jj,jk) 
     202                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk) 
    284203                  END DO 
    285204               END DO 
    286205            END DO 
    287206         ENDIF 
    288          pht_adv(:) = ptr_vj( zt2(:,:,:) ) 
    289          pst_adv(:) = ptr_vj( zs2(:,:,:) ) 
    290       ENDIF 
     207         IF( ktra == jp_tem)   pht_adv(:) = ptr_vj( zwy(:,:,:) ) 
     208         IF( ktra == jp_sal)   pst_adv(:) = ptr_vj( zwy(:,:,:) ) 
     209      ENDIF 
     210 
    291211 
    292212      ! II. Vertical advective fluxes 
     
    296216      ! interior values 
    297217      DO jk = 2, jpkm1 
    298          zt1(:,:,jk) = tmask(:,:,jk) * ( tb(:,:,jk-1) - tb(:,:,jk) ) 
    299          zs1(:,:,jk) = tmask(:,:,jk) * ( sb(:,:,jk-1) - sb(:,:,jk) ) 
     218         zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1) - ptb(:,:,jk) ) 
    300219      END DO 
    301220      ! surface & bottom boundary conditions 
    302       zt1 (:,:, 1 ) = 0.e0    ;    zt1 (:,:,jpk) = 0.e0 
    303       zs1 (:,:, 1 ) = 0.e0    ;    zs1 (:,:,jpk) = 0.e0 
     221      zwx (:,:, 1 ) = 0.e0    ;    zwx (:,:,jpk) = 0.e0 
    304222 
    305223      ! Slopes 
     
    307225         DO jj = 1, jpj 
    308226            DO ji = 1, jpi 
    309                ztp1(ji,jj,jk) =                    ( zt1(ji,jj,jk) + zt1(ji,jj,jk+1) )   & 
    310                   &           * ( 0.25 + SIGN( 0.25, zt1(ji,jj,jk) * zt1(ji,jj,jk+1) ) ) 
    311                zsp1(ji,jj,jk) =                    ( zs1(ji,jj,jk) + zs1(ji,jj,jk+1) )   & 
    312                   &           * ( 0.25 + SIGN( 0.25, zs1(ji,jj,jk) * zs1(ji,jj,jk+1) ) ) 
     227               zslpx(ji,jj,jk) =                    ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) )   & 
     228                  &            * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) ) 
    313229            END DO 
    314230         END DO 
     
    316232 
    317233      ! Slopes limitation 
    318       ! interior values 
    319       DO jk = 2, jpkm1 
     234      DO jk = 2, jpkm1        ! interior values 
    320235         DO jj = 1, jpj 
    321236            DO ji = 1, jpi 
    322                ztp1(ji,jj,jk) = SIGN( 1., ztp1(ji,jj,jk) )   & 
    323                   &           * MIN(    ABS( ztp1(ji,jj,jk  ) ),   & 
    324                   &                  2.*ABS( zt1 (ji,jj,jk+1) ),   & 
    325                   &                  2.*ABS( zt1 (ji,jj,jk  ) ) ) 
    326                zsp1(ji,jj,jk) = SIGN( 1., zsp1(ji,jj,jk) )   & 
    327                   &           * MIN(    ABS( zsp1(ji,jj,jk  ) ),   & 
    328                   &                  2.*ABS( zs1 (ji,jj,jk+1) ),   & 
    329                   &                  2.*ABS( zs1 (ji,jj,jk  ) ) ) 
    330             END DO 
    331          END DO 
    332       END DO 
    333       ! surface values 
    334       ztp1(:,:,1) = 0.e0 
    335       zsp1(:,:,1) = 0.e0 
     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 
    336245 
    337246      ! vertical advective flux 
    338       ! interior values 
    339       DO jk = 1, jpkm1 
     247      DO jk = 1, jpkm1        ! interior values 
    340248         zstep  = z2 * rdttra(jk) 
    341249         DO jj = 2, jpjm1       
     
    344252               z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 
    345253               zalpha = 0.5 + z0w 
    346                zw  = z0w - 0.5 * pwn(ji,jj,jk+1)*zstep / fse3w(ji,jj,jk+1) 
    347                zzt1 = tb(ji,jj,jk+1) + zw*ztp1(ji,jj,jk+1) 
    348                zzt2 = tb(ji,jj,jk  ) + zw*ztp1(ji,jj,jk  ) 
    349                zzs1 = sb(ji,jj,jk+1) + zw*zsp1(ji,jj,jk+1) 
    350                zzs2 = sb(ji,jj,jk  ) + zw*zsp1(ji,jj,jk  ) 
    351                zt1(ji,jj,jk+1) = zew * ( zalpha * zzt1 + (1.-zalpha)*zzt2 ) 
    352                zs1(ji,jj,jk+1) = zew * ( zalpha * zzs1 + (1.-zalpha)*zzs2 ) 
    353             END DO 
    354          END DO 
    355       END DO 
    356       ! surface values 
     254               zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1) 
     255               zzwx = ptb(ji,jj,jk+1) + zw * zslpx(ji,jj,jk+1) 
     256               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 
    357262      IF( lk_dynspg_rl .OR. lk_vvl) THEN                ! rigid lid or variable volume: flux set to zero 
    358          zt1(:,:, 1 ) = 0.e0 
    359          zs1(:,:, 1 ) = 0.e0 
     263         zwx(:,:, 1 ) = 0.e0                                 ! surface  
     264         zwx(:,:,jpk) = 0.e0                                 ! bottom  
    360265      ELSE                                              ! free surface 
    361          zt1(:,:, 1 ) = pwn(:,:,1) * tb(:,:,1) 
    362          zs1(:,:, 1 ) = pwn(:,:,1) * sb(:,:,1) 
    363       ENDIF 
    364  
    365       ! bottom values 
    366       zt1(:,:,jpk) = 0.e0 
    367       zs1(:,:,jpk) = 0.e0 
    368  
     266         zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1)              ! Surface 
     267         zwx(:,:,jpk) = 0.e0                                 ! bottom  
     268 
     269      ENDIF 
    369270 
    370271      ! Compute & add the vertical advective trend 
    371  
    372272      DO jk = 1, jpkm1 
    373273         DO jj = 2, jpjm1       
     
    375275               zbtr = 1. / fse3t(ji,jj,jk) 
    376276               ! horizontal advective trends 
    377                zta = - zbtr * ( zt1(ji,jj,jk) - zt1(ji,jj,jk+1) ) 
    378                zsa = - zbtr * ( zs1(ji,jj,jk) - zs1(ji,jj,jk+1) ) 
     277               zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 
    379278               ! add it to the general tracer trends 
    380                ta(ji,jj,jk) =  ta(ji,jj,jk) + zta 
    381                sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa 
     279               pta(ji,jj,jk) =  pta(ji,jj,jk) + zta 
    382280            END DO 
    383281         END DO 
     
    386284      ! Save the vertical advective trends for diagnostic 
    387285      ! ------------------------------------------------- 
    388       IF( l_trdtra )   THEN 
    389          ! Recompute the vertical advection zta & zsa trends computed  
    390          ! at the step 2. above in making the difference between the new  
    391          ! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract 
    392          ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends 
    393  
    394          DO jk = 1, jpkm1 
    395             DO jj = 2, jpjm1 
    396                DO ji = fs_2, fs_jpim1   ! vector opt. 
    397 #if defined key_zco 
    398                   zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) ) 
    399                   z_hdivn_x = e2u(ji,jj)*pun(ji,jj,jk) - e2u(ji-1,jj)*pun(ji-1,jj,jk) 
    400                   z_hdivn_y = e1v(ji,jj)*pvn(ji,jj,jk) - e1v(ji,jj-1)*pvn(ji,jj-1,jk) 
    401 #else 
    402                   zbtr      = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    403                   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) 
    404                   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) 
    405 #endif 
    406                   z_hdivn   = (z_hdivn_x + z_hdivn_y) * zbtr 
    407                   ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn  
    408                   ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn 
    409                END DO 
    410             END DO 
    411          END DO 
    412          CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt) 
    413          ! 
    414       ENDIF 
    415  
    416       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' muscl zad  - Ta: ', mask1=tmask ,   & 
    417          &                       tab3d_2=sa, clinfo2=             ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     286      IF( l_trdtra )   CALL trd_tra_adv( kt, ktra, jpt_trd_zad, cdtype, zwx, pwn, ptb ) 
     287 
     288      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pta, clinfo1=' muscl - zad: ', mask1=tmask, clinfo3=cdtype ) 
    418289      ! 
    419290   END SUBROUTINE tra_adv_muscl 
Note: See TracChangeset for help on using the changeset viewer.