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 12377 for NEMO/trunk/src/OCE/TRA/traadv_fct.F90 – NEMO

Ignore:
Timestamp:
2020-02-12T15:39:06+01:00 (4 years ago)
Author:
acc
Message:

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 

  • NEMO/trunk/src/OCE/TRA/traadv_fct.F90

    r12055 r12377  
    4545 
    4646   !! * Substitutions 
    47 #  include "vectopt_loop_substitute.h90" 
     47#  include "do_loop_substitute.h90" 
    4848   !!---------------------------------------------------------------------- 
    4949   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    5353CONTAINS 
    5454 
    55    SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn,       & 
    56       &                                              ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v ) 
     55   SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pU, pV, pW,       & 
     56      &                    Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) 
    5757      !!---------------------------------------------------------------------- 
    5858      !!                  ***  ROUTINE tra_adv_fct  *** 
    … …  
    6666      !!               - corrected flux (monotonic correction)  
    6767      !! 
    68       !! ** Action : - update pta  with the now advective tracer trends 
     68      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends 
    6969      !!             - send trends to trdtra module for further diagnostics (l_trdtra=T) 
    70       !!             - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T) 
    71       !!---------------------------------------------------------------------- 
    72       INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index 
    73       INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index 
    74       CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
    75       INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers 
    76       INTEGER                              , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4) 
    77       INTEGER                              , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4) 
    78       REAL(wp)                             , INTENT(in   ) ::   p2dt            ! tracer time-step 
    79       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components 
    80       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields 
    81       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
     70      !!             - poleward advective heat and salt transport (ln_diaptr=T) 
     71      !!---------------------------------------------------------------------- 
     72      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index 
     73      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     74      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index 
     75      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
     76      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers 
     77      INTEGER                                  , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4) 
     78      INTEGER                                  , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4) 
     79      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
     80      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components 
     81      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    8282      ! 
    8383      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices   
    … …  
    101101      l_ptr = .FALSE. 
    102102      ll_zAimp = .FALSE. 
    103       IF( ( cdtype =='TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )       l_trd = .TRUE. 
    104       IF(   cdtype =='TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE.  
    105       IF(   cdtype =='TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.  & 
     103      IF( ( cdtype == 'TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
     104      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    l_ptr = .TRUE.  
     105      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.  & 
    106106         &                         iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE. 
    107107      ! 
    … …  
    128128      IF( ll_zAimp ) THEN 
    129129         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 
    130          DO jk = 1, jpkm1 
    131             DO jj = 2, jpjm1 
    132                DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.) 
    133                   zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t_a(ji,jj,jk) 
    134                   zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t_a(ji,jj,jk) 
    135                   zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t_a(ji,jj,jk) 
    136                END DO 
    137             END DO 
    138          END DO 
     130         DO_3D_00_00( 1, jpkm1 ) 
     131            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t(ji,jj,jk,Krhs) 
     132            zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t(ji,jj,jk,Krhs) 
     133            zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) 
     134         END_3D 
    139135      END IF 
    140136      ! 
    … …  
    143139         !        !==  upstream advection with initial mass fluxes & intermediate update  ==! 
    144140         !                    !* upstream tracer flux in the i and j direction  
    145          DO jk = 1, jpkm1 
    146             DO jj = 1, jpjm1 
    147                DO ji = 1, fs_jpim1   ! vector opt. 
    148                   ! upstream scheme 
    149                   zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 
    150                   zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 
    151                   zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 
    152                   zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 
    153                   zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn) ) 
    154                   zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn) ) 
    155                END DO 
    156             END DO 
    157          END DO 
     141         DO_3D_10_10( 1, jpkm1 ) 
     142            ! upstream scheme 
     143            zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) 
     144            zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) ) 
     145            zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) ) 
     146            zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) ) 
     147            zwx(ji,jj,jk) = 0.5 * ( zfp_ui * pt(ji,jj,jk,jn,Kbb) + zfm_ui * pt(ji+1,jj  ,jk,jn,Kbb) ) 
     148            zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) ) 
     149         END_3D 
    158150         !                    !* upstream tracer flux in the k direction *! 
    159          DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
    160             DO jj = 1, jpj 
    161                DO ji = 1, jpi 
    162                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    163                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    164                   zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
    165                END DO 
    166             END DO 
    167          END DO 
     151         DO_3D_11_11( 2, jpkm1 ) 
     152            zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 
     153            zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 
     154            zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk) 
     155         END_3D 
    168156         IF( ln_linssh ) THEN    ! top ocean value (only in linear free surface as zwz has been w-masked) 
    169157            IF( ln_isfcav ) THEN             ! top of the ice-shelf cavities and at the ocean surface 
    170                DO jj = 1, jpj 
    171                   DO ji = 1, jpi 
    172                      zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    173                   END DO 
    174                END DO    
     158               DO_2D_11_11 
     159                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface  
     160               END_2D 
    175161            ELSE                             ! no cavities: only at the ocean surface 
    176                zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 
     162               zwz(:,:,1) = pW(:,:,1) * pt(:,:,1,jn,Kbb) 
    177163            ENDIF 
    178164         ENDIF 
    179165         !                
    180          DO jk = 1, jpkm1     !* trend and after field with monotonic scheme 
    181             DO jj = 2, jpjm1 
    182                DO ji = fs_2, fs_jpim1   ! vector opt. 
    183                   !                             ! total intermediate advective trends 
    184                   ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    185                      &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    186                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
    187                   !                             ! update and guess with monotonic sheme 
    188                   pta(ji,jj,jk,jn) =                     pta(ji,jj,jk,jn) +        ztra   / e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
    189                   zwi(ji,jj,jk)    = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    190                END DO 
    191             END DO 
    192          END DO 
     166         DO_3D_00_00( 1, jpkm1 ) 
     167            !                             ! total intermediate advective trends 
     168            ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     169               &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     170               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     171            !                             ! update and guess with monotonic sheme 
     172            pt(ji,jj,jk,jn,Krhs) =                     pt(ji,jj,jk,jn,Krhs) +        ztra   / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) 
     173            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
     174         END_3D 
    193175          
    194176         IF ( ll_zAimp ) THEN 
    … …  
    196178            ! 
    197179            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 
    198             DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
    199                DO jj = 2, jpjm1 
    200                   DO ji = fs_2, fs_jpim1   ! vector opt.   
    201                      zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
    202                      zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
    203                      ztw(ji,jj,jk) =  0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
    204                      zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 
    205                   END DO 
    206                END DO 
    207             END DO 
    208             DO jk = 1, jpkm1 
    209                DO jj = 2, jpjm1 
    210                   DO ji = fs_2, fs_jpim1   ! vector opt.   
    211                      pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
    212                         &                                  * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    213                   END DO 
    214                END DO 
    215             END DO 
     180            DO_3D_00_00( 2, jpkm1 ) 
     181               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     182               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     183               ztw(ji,jj,jk) =  0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     184               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 
     185            END_3D 
     186            DO_3D_00_00( 1, jpkm1 ) 
     187               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
     188                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     189            END_3D 
    216190            ! 
    217191         END IF 
    … …  
    228202         ! 
    229203         CASE(  2  )                   !- 2nd order centered 
    230             DO jk = 1, jpkm1 
    231                DO jj = 1, jpjm1 
    232                   DO ji = 1, fs_jpim1   ! vector opt. 
    233                      zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk) 
    234                      zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk) 
    235                   END DO 
    236                END DO 
    237             END DO 
     204            DO_3D_10_10( 1, jpkm1 ) 
     205               zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk) 
     206               zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk) 
     207            END_3D 
    238208            ! 
    239209         CASE(  4  )                   !- 4th order centered 
    … …  
    241211            zltv(:,:,jpk) = 0._wp 
    242212            DO jk = 1, jpkm1                 ! Laplacian 
    243                DO jj = 1, jpjm1                    ! 1st derivative (gradient) 
    244                   DO ji = 1, fs_jpim1   ! vector opt. 
    245                      ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    246                      ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    247                   END DO 
    248                END DO 
    249                DO jj = 2, jpjm1                    ! 2nd derivative * 1/ 6 
    250                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    251                      zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6 
    252                      zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6 
    253                   END DO 
    254                END DO 
     213               DO_2D_10_10 
     214                  ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 
     215                  ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
     216               END_2D 
     217               DO_2D_00_00 
     218                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6 
     219                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6 
     220               END_2D 
    255221            END DO 
    256222            CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1. , zltv, 'T', 1. )   ! Lateral boundary cond. (unchanged sgn) 
    257223            ! 
    258             DO jk = 1, jpkm1                 ! Horizontal advective fluxes 
    259                DO jj = 1, jpjm1 
    260                   DO ji = 1, fs_jpim1   ! vector opt. 
    261                      zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points 
    262                      zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) 
    263                      !                                                  ! C4 minus upstream advective fluxes  
    264                      zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 
    265                      zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) 
    266                   END DO 
    267                END DO 
    268             END DO          
     224            DO_3D_10_10( 1, jpkm1 ) 
     225               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points 
     226               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) 
     227               !                                                  ! C4 minus upstream advective fluxes  
     228               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 
     229               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) 
     230            END_3D 
    269231            ! 
    270232         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested 
    271233            ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero 
    272234            ztv(:,:,jpk) = 0._wp 
    273             DO jk = 1, jpkm1                 ! 1st derivative (gradient) 
    274                DO jj = 1, jpjm1 
    275                   DO ji = 1, fs_jpim1   ! vector opt. 
    276                      ztu(ji,jj,jk) = ( ptn(ji+1,jj  ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    277                      ztv(ji,jj,jk) = ( ptn(ji  ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    278                   END DO 
    279                END DO 
    280             END DO 
     235            DO_3D_10_10( 1, jpkm1 ) 
     236               ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 
     237               ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
     238            END_3D 
    281239            CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1. , ztv, 'V', -1. )   ! Lateral boundary cond. (unchanged sgn) 
    282240            ! 
    283             DO jk = 1, jpkm1                 ! Horizontal advective fluxes 
    284                DO jj = 2, jpjm1 
    285                   DO ji = 2, fs_jpim1   ! vector opt. 
    286                      zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn)   ! 2 x C2 interpolation of T at u- & v-points (x2) 
    287                      zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) 
    288                      !                                                  ! C4 interpolation of T at u- & v-points (x2) 
    289                      zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) ) 
    290                      zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) ) 
    291                      !                                                  ! C4 minus upstream advective fluxes  
    292                      zwx(ji,jj,jk) =  0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 
    293                      zwy(ji,jj,jk) =  0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 
    294                   END DO 
    295                END DO 
    296             END DO 
     241            DO_3D_00_00( 1, jpkm1 ) 
     242               zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj  ,jk,jn,Kmm)   ! 2 x C2 interpolation of T at u- & v-points (x2) 
     243               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) 
     244               !                                                  ! C4 interpolation of T at u- & v-points (x2) 
     245               zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) ) 
     246               zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) ) 
     247               !                                                  ! C4 minus upstream advective fluxes  
     248               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 
     249               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 
     250            END_3D 
    297251            ! 
    298252         END SELECT 
    … …  
    301255         ! 
    302256         CASE(  2  )                   !- 2nd order centered 
    303             DO jk = 2, jpkm1     
    304                DO jj = 2, jpjm1 
    305                   DO ji = fs_2, fs_jpim1 
    306                      zwz(ji,jj,jk) =  (  pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )   & 
    307                         &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk) 
    308                   END DO 
    309                END DO 
    310             END DO 
     257            DO_3D_00_00( 2, jpkm1 ) 
     258               zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   & 
     259                  &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk) 
     260            END_3D 
    311261            ! 
    312262         CASE(  4  )                   !- 4th order COMPACT 
    313             CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw )   ! zwt = COMPACT interpolation of T at w-point 
    314             DO jk = 2, jpkm1 
    315                DO jj = 2, jpjm1 
    316                   DO ji = fs_2, fs_jpim1 
    317                      zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
    318                   END DO 
    319                END DO 
    320             END DO 
     263            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )   ! zwt = COMPACT interpolation of T at w-point 
     264            DO_3D_00_00( 2, jpkm1 ) 
     265               zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
     266            END_3D 
    321267            ! 
    322268         END SELECT 
    … …  
    326272         !          
    327273         IF ( ll_zAimp ) THEN 
    328             DO jk = 1, jpkm1     !* trend and after field with monotonic scheme 
    329                DO jj = 2, jpjm1 
    330                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    331                      !                             ! total intermediate advective trends 
    332                      ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    333                         &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    334                         &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
    335                      ztw(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    336                   END DO 
    337                END DO 
    338             END DO 
     274            DO_3D_00_00( 1, jpkm1 ) 
     275               !                             ! total intermediate advective trends 
     276               ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     277                  &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     278                  &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     279               ztw(ji,jj,jk)  = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
     280            END_3D 
    339281            ! 
    340282            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 
    341283            ! 
    342             DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
    343                DO jj = 2, jpjm1 
    344                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    345                      zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
    346                      zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
    347                      zwz(ji,jj,jk) =  zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
    348                   END DO 
    349                END DO 
    350             END DO 
     284            DO_3D_00_00( 2, jpkm1 ) 
     285               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     286               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     287               zwz(ji,jj,jk) =  zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     288            END_3D 
    351289         END IF 
    352290         ! 
    … …  
    355293         !        !==  monotonicity algorithm  ==! 
    356294         ! 
    357          CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 
     295         CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx, zwy, zwz, zwi, p2dt ) 
    358296         ! 
    359297         !        !==  final trend with corrected fluxes  ==! 
    360298         ! 
    361          DO jk = 1, jpkm1 
    362             DO jj = 2, jpjm1 
    363                DO ji = fs_2, fs_jpim1   ! vector opt.   
    364                   ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    365                      &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    366                      &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
    367                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) 
    368                   zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 
    369                END DO 
    370             END DO 
    371          END DO 
     299         DO_3D_00_00( 1, jpkm1 ) 
     300            ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     301               &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     302               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     303            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) 
     304            zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
     305         END_3D 
    372306         ! 
    373307         IF ( ll_zAimp ) THEN 
    374308            ! 
    375309            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 
    376             DO jk = 2, jpkm1        ! Interior value ( multiplied by wmask) 
    377                DO jj = 2, jpjm1 
    378                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    379                      zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
    380                      zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
    381                      ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
    382                      zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 
    383                   END DO 
    384                END DO 
    385             END DO 
    386             DO jk = 1, jpkm1 
    387                DO jj = 2, jpjm1 
    388                   DO ji = fs_2, fs_jpim1   ! vector opt.   
    389                      pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
    390                         &                                * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 
    391                   END DO 
    392                END DO 
    393             END DO 
     310            DO_3D_00_00( 2, jpkm1 ) 
     311               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     312               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     313               ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 
     314               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 
     315            END_3D 
     316            DO_3D_00_00( 1, jpkm1 ) 
     317               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
     318                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     319            END_3D 
    394320         END IF          
    395321         ! 
    … …  
    400326            ! 
    401327            IF( l_trd ) THEN              ! trend diagnostics 
    402                CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 
    403                CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 
    404                CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 
     328               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) ) 
     329               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) ) 
     330               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) ) 
    405331            ENDIF 
    406332            !                             ! heat/salt transport 
    … …  
    428354 
    429355 
    430    SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt ) 
     356   SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt ) 
    431357      !!--------------------------------------------------------------------- 
    432358      !!                    ***  ROUTINE nonosc  *** 
    … …  
    441367      !!       in-space based differencing for fluid 
    442368      !!---------------------------------------------------------------------- 
     369      INTEGER                          , INTENT(in   ) ::   Kmm             ! time level index  
    443370      REAL(wp)                         , INTENT(in   ) ::   p2dt            ! tracer time-step 
    444371      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pbef, paft      ! before & after field 
    … …  
    466393      DO jk = 1, jpkm1 
    467394         ikm1 = MAX(jk-1,1) 
    468          DO jj = 2, jpjm1 
    469             DO ji = fs_2, fs_jpim1   ! vector opt. 
    470  
    471                ! search maximum in neighbourhood 
    472                zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 
    473                   &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   & 
    474                   &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   & 
    475                   &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  ) 
    476  
    477                ! search minimum in neighbourhood 
    478                zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   & 
    479                   &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   & 
    480                   &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   & 
    481                   &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  ) 
    482  
    483                ! positive part of the flux 
    484                zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   & 
    485                   & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   & 
    486                   & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) ) 
    487  
    488                ! negative part of the flux 
    489                zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   & 
    490                   & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   & 
    491                   & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
    492  
    493                ! up & down beta terms 
    494                zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 
    495                zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 
    496                zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt 
    497             END DO 
    498          END DO 
     395         DO_2D_00_00 
     396 
     397            ! search maximum in neighbourhood 
     398            zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 
     399               &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   & 
     400               &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   & 
     401               &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  ) 
     402 
     403            ! search minimum in neighbourhood 
     404            zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   & 
     405               &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   & 
     406               &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   & 
     407               &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  ) 
     408 
     409            ! positive part of the flux 
     410            zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   & 
     411               & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   & 
     412               & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) ) 
     413 
     414            ! negative part of the flux 
     415            zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   & 
     416               & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   & 
     417               & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
     418 
     419            ! up & down beta terms 
     420            zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt 
     421            zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 
     422            zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt 
     423         END_2D 
    499424      END DO 
    500425      CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1. , zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
    … …  
    502427      ! 3. monotonic flux in the i & j direction (paa & pbb) 
    503428      ! ---------------------------------------- 
    504       DO jk = 1, jpkm1 
    505          DO jj = 2, jpjm1 
    506             DO ji = fs_2, fs_jpim1   ! vector opt. 
    507                zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
    508                zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
    509                zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) ) 
    510                paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
    511  
    512                zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
    513                zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
    514                zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 
    515                pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
    516  
    517       ! monotonic flux in the k direction, i.e. pcc 
    518       ! ------------------------------------------- 
    519                za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 
    520                zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 
    521                zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) ) 
    522                pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 
    523             END DO 
    524          END DO 
    525       END DO 
     429      DO_3D_00_00( 1, jpkm1 ) 
     430         zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
     431         zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
     432         zcu =       ( 0.5  + SIGN( 0.5 , paa(ji,jj,jk) ) ) 
     433         paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
     434 
     435         zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
     436         zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
     437         zcv =       ( 0.5  + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 
     438         pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
     439 
     440! monotonic flux in the k direction, i.e. pcc 
     441! ------------------------------------------- 
     442         za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 
     443         zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 
     444         zc =       ( 0.5  + SIGN( 0.5 , pcc(ji,jj,jk+1) ) ) 
     445         pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 
     446      END_3D 
    526447      CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1. , pbb, 'V', -1. )   ! lateral boundary condition (changed sign) 
    527448      ! 
    … …  
    544465      !!---------------------------------------------------------------------- 
    545466       
    546       DO jk = 3, jpkm1        !==  build the three diagonal matrix  ==! 
    547          DO jj = 1, jpj 
    548             DO ji = 1, jpi 
    549                zwd (ji,jj,jk) = 4._wp 
    550                zwi (ji,jj,jk) = 1._wp 
    551                zws (ji,jj,jk) = 1._wp 
    552                zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
    553                ! 
    554                IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom 
    555                   zwd (ji,jj,jk) = 1._wp 
    556                   zwi (ji,jj,jk) = 0._wp 
    557                   zws (ji,jj,jk) = 0._wp 
    558                   zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )     
    559                ENDIF 
    560             END DO 
    561          END DO 
    562       END DO 
    563       ! 
    564       jk = 2                                          ! Switch to second order centered at top 
    565       DO jj = 1, jpj 
    566          DO ji = 1, jpi 
     467      DO_3D_11_11( 3, jpkm1 ) 
     468         zwd (ji,jj,jk) = 4._wp 
     469         zwi (ji,jj,jk) = 1._wp 
     470         zws (ji,jj,jk) = 1._wp 
     471         zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
     472         ! 
     473         IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom 
    567474            zwd (ji,jj,jk) = 1._wp 
    568475            zwi (ji,jj,jk) = 0._wp 
    569476            zws (ji,jj,jk) = 0._wp 
    570             zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
    571          END DO 
    572       END DO    
     477            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )     
     478         ENDIF 
     479      END_3D 
     480      ! 
     481      jk = 2                                          ! Switch to second order centered at top 
     482      DO_2D_11_11 
     483         zwd (ji,jj,jk) = 1._wp 
     484         zwi (ji,jj,jk) = 0._wp 
     485         zws (ji,jj,jk) = 0._wp 
     486         zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
     487      END_2D 
    573488      ! 
    574489      !                       !==  tridiagonal solve  ==! 
    575       DO jj = 1, jpj                ! first recurrence 
    576          DO ji = 1, jpi 
    577             zwt(ji,jj,2) = zwd(ji,jj,2) 
    578          END DO 
    579       END DO 
    580       DO jk = 3, jpkm1 
    581          DO jj = 1, jpj 
    582             DO ji = 1, jpi 
    583                zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
    584             END DO 
    585          END DO 
    586       END DO 
    587       ! 
    588       DO jj = 1, jpj                ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    589          DO ji = 1, jpi 
    590             pt_out(ji,jj,2) = zwrm(ji,jj,2) 
    591          END DO 
    592       END DO 
    593       DO jk = 3, jpkm1 
    594          DO jj = 1, jpj 
    595             DO ji = 1, jpi 
    596                pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
    597             END DO 
    598          END DO 
    599       END DO 
    600  
    601       DO jj = 1, jpj                ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 
    602          DO ji = 1, jpi 
    603             pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
    604          END DO 
    605       END DO 
    606       DO jk = jpk-2, 2, -1 
    607          DO jj = 1, jpj 
    608             DO ji = 1, jpi 
    609                pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    610             END DO 
    611          END DO 
    612       END DO 
     490      DO_2D_11_11 
     491         zwt(ji,jj,2) = zwd(ji,jj,2) 
     492      END_2D 
     493      DO_3D_11_11( 3, jpkm1 ) 
     494         zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
     495      END_3D 
     496      ! 
     497      DO_2D_11_11 
     498         pt_out(ji,jj,2) = zwrm(ji,jj,2) 
     499      END_2D 
     500      DO_3D_11_11( 3, jpkm1 ) 
     501         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
     502      END_3D 
     503 
     504      DO_2D_11_11 
     505         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
     506      END_2D 
     507      DO_3DS_11_11( jpk-2, 2, -1 ) 
     508         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
     509      END_3D 
    613510      !     
    614511   END SUBROUTINE interp_4th_cpt_org 
    … …  
    633530      !                      !==  build the three diagonal matrix & the RHS  ==! 
    634531      ! 
    635       DO jk = 3, jpkm1                 ! interior (from jk=3 to jpk-1) 
    636          DO jj = 2, jpjm1 
    637             DO ji = fs_2, fs_jpim1 
    638                zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal 
    639                zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal 
    640                zws (ji,jj,jk) =         wmask(ji,jj,jk)                         ! upper diagonal 
    641                zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk)                     &   ! RHS 
    642                   &           *       ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 
    643             END DO 
    644          END DO 
    645       END DO 
     532      DO_3D_00_00( 3, jpkm1 ) 
     533         zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal 
     534         zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal 
     535         zws (ji,jj,jk) =         wmask(ji,jj,jk)                         ! upper diagonal 
     536         zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk)                     &   ! RHS 
     537            &           *       ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 
     538      END_3D 
    646539      ! 
    647540!!gm 
    … …  
    656549      END IF 
    657550      ! 
    658       DO jj = 2, jpjm1                 ! 2nd order centered at top & bottom 
    659          DO ji = fs_2, fs_jpim1 
    660             ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point 
    661             ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point 
    662             ! 
    663             zwd (ji,jj,ikt) = 1._wp          ! top 
    664             zwi (ji,jj,ikt) = 0._wp 
    665             zws (ji,jj,ikt) = 0._wp 
    666             zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) ) 
    667             ! 
    668             zwd (ji,jj,ikb) = 1._wp          ! bottom 
    669             zwi (ji,jj,ikb) = 0._wp 
    670             zws (ji,jj,ikb) = 0._wp 
    671             zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )             
    672          END DO 
    673       END DO    
     551      DO_2D_00_00 
     552         ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point 
     553         ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point 
     554         ! 
     555         zwd (ji,jj,ikt) = 1._wp          ! top 
     556         zwi (ji,jj,ikt) = 0._wp 
     557         zws (ji,jj,ikt) = 0._wp 
     558         zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) ) 
     559         ! 
     560         zwd (ji,jj,ikb) = 1._wp          ! bottom 
     561         zwi (ji,jj,ikb) = 0._wp 
     562         zws (ji,jj,ikb) = 0._wp 
     563         zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) )             
     564      END_2D 
    674565      ! 
    675566      !                       !==  tridiagonal solver  ==! 
    676567      ! 
    677       DO jj = 2, jpjm1              !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1 
    678          DO ji = fs_2, fs_jpim1 
    679             zwt(ji,jj,2) = zwd(ji,jj,2) 
    680          END DO 
    681       END DO 
    682       DO jk = 3, jpkm1 
    683          DO jj = 2, jpjm1 
    684             DO ji = fs_2, fs_jpim1 
    685                zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
    686             END DO 
    687          END DO 
    688       END DO 
    689       ! 
    690       DO jj = 2, jpjm1              !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    691          DO ji = fs_2, fs_jpim1 
    692             pt_out(ji,jj,2) = zwrm(ji,jj,2) 
    693          END DO 
    694       END DO 
    695       DO jk = 3, jpkm1 
    696          DO jj = 2, jpjm1 
    697             DO ji = fs_2, fs_jpim1 
    698                pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
    699             END DO 
    700          END DO 
    701       END DO 
    702  
    703       DO jj = 2, jpjm1              !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk 
    704          DO ji = fs_2, fs_jpim1 
    705             pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
    706          END DO 
    707       END DO 
    708       DO jk = jpk-2, 2, -1 
    709          DO jj = 2, jpjm1 
    710             DO ji = fs_2, fs_jpim1 
    711                pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    712             END DO 
    713          END DO 
    714       END DO 
     568      DO_2D_00_00 
     569         zwt(ji,jj,2) = zwd(ji,jj,2) 
     570      END_2D 
     571      DO_3D_00_00( 3, jpkm1 ) 
     572         zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
     573      END_3D 
     574      ! 
     575      DO_2D_00_00 
     576         pt_out(ji,jj,2) = zwrm(ji,jj,2) 
     577      END_2D 
     578      DO_3D_00_00( 3, jpkm1 ) 
     579         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
     580      END_3D 
     581 
     582      DO_2D_00_00 
     583         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
     584      END_2D 
     585      DO_3DS_00_00( jpk-2, 2, -1 ) 
     586         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
     587      END_3D 
    715588      !     
    716589   END SUBROUTINE interp_4th_cpt 
    … …  
    749622      kstart =  1  + klev 
    750623      ! 
    751       DO jj = 2, jpjm1              !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1 
    752          DO ji = fs_2, fs_jpim1 
    753             zwt(ji,jj,kstart) = pD(ji,jj,kstart) 
    754          END DO 
    755       END DO 
    756       DO jk = kstart+1, jpkm1 
    757          DO jj = 2, jpjm1 
    758             DO ji = fs_2, fs_jpim1 
    759                zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
    760             END DO 
    761          END DO 
    762       END DO 
    763       ! 
    764       DO jj = 2, jpjm1              !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    765          DO ji = fs_2, fs_jpim1 
    766             pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 
    767          END DO 
    768       END DO 
    769       DO jk = kstart+1, jpkm1 
    770          DO jj = 2, jpjm1 
    771             DO ji = fs_2, fs_jpim1 
    772                pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
    773             END DO 
    774          END DO 
    775       END DO 
    776  
    777       DO jj = 2, jpjm1              !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk 
    778          DO ji = fs_2, fs_jpim1 
    779             pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
    780          END DO 
    781       END DO 
    782       DO jk = jpk-2, kstart, -1 
    783          DO jj = 2, jpjm1 
    784             DO ji = fs_2, fs_jpim1 
    785                pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    786             END DO 
    787          END DO 
    788       END DO 
     624      DO_2D_00_00 
     625         zwt(ji,jj,kstart) = pD(ji,jj,kstart) 
     626      END_2D 
     627      DO_3D_00_00( kstart+1, jpkm1 ) 
     628         zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
     629      END_3D 
     630      ! 
     631      DO_2D_00_00 
     632         pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 
     633      END_2D 
     634      DO_3D_00_00( kstart+1, jpkm1 ) 
     635         pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)              
     636      END_3D 
     637 
     638      DO_2D_00_00 
     639         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
     640      END_2D 
     641      DO_3DS_00_00( jpk-2, kstart, -1 ) 
     642         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
     643      END_3D 
    789644      ! 
    790645   END SUBROUTINE tridia_solver 
Note: See TracChangeset for help on using the changeset viewer.