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 14986 for NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA – NEMO

Ignore:
Timestamp:
2021-06-14T13:34:08+02:00 (3 years ago)
Author:
sparonuz
Message:

Merge trunk -r14984:HEAD

Location:
NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA
Files:
3 added
2 deleted
21 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/eosbn2.F90

    r14219 r14986  
    578578 
    579579   SUBROUTINE eos_insitu_pot_2d( pts, prhop ) 
     580      !! 
     581      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
     582      !                                                     ! 2 : salinity               [psu] 
     583      REAL(dp), DIMENSION(:,:)  , INTENT(  out) ::   prhop  ! potential density (surface referenced) 
     584      !! 
     585      CALL eos_insitu_pot_2d_t( pts, is_tile(pts), prhop, is_tile(prhop) ) 
     586   END SUBROUTINE eos_insitu_pot_2d 
     587 
     588 
     589   SUBROUTINE eos_insitu_pot_2d_t( pts, ktts, prhop, ktrhop ) 
    580590      !!---------------------------------------------------------------------- 
    581591      !!                  ***  ROUTINE eos_insitu_pot  *** 
    … …  
    590600      !! 
    591601      !!---------------------------------------------------------------------- 
    592       REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
     602      INTEGER                              , INTENT(in   ) ::   ktts, ktrhop 
     603      REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius] 
    593604      !                                                                ! 2 : salinity               [psu] 
    594       REAL(dp), DIMENSION(jpi,jpj     ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
     605      REAL(dp), DIMENSION(A2D_T(ktrhop)   ), INTENT(  out) ::   prhop  ! potential density (surface referenced) 
    595606      ! 
    596607      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
    … …  
    607618      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    608619         ! 
    609             DO_2D( 1, 1, 1, 1 ) 
    610                ! 
    611                zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
    612                zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    613                ztm = tmask(ji,jj,1)                                         ! tmask 
    614                ! 
    615                zn0 = (((((EOS060*zt   & 
    616                   &   + EOS150*zs+EOS050)*zt   & 
    617                   &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    618                   &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    619                   &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    620                   &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    621                   &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    622                   ! 
    623                ! 
    624                prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface 
    625                ! 
    626             END_2D 
     620         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     621            ! 
     622            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature 
     623            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     624            ztm = tmask(ji,jj,1)                                         ! tmask 
     625            ! 
     626            zn0 = (((((EOS060*zt   & 
     627               &   + EOS150*zs+EOS050)*zt   & 
     628               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     629               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     630               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     631               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     632               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     633               ! 
     634            ! 
     635            prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface 
     636            ! 
     637         END_2D 
    627638 
    628639      CASE( np_seos )                !==  simplified EOS  ==! 
    629640         ! 
    630          DO_2D( 1, 1, 1, 1 ) 
     641         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    631642            zt  = pts  (ji,jj,jp_tem) - 10._wp 
    632643            zs  = pts  (ji,jj,jp_sal) - 35._wp 
    … …  
    647658      IF( ln_timing )   CALL timing_stop('eos-pot') 
    648659      ! 
    649    END SUBROUTINE eos_insitu_pot_2d 
     660   END SUBROUTINE eos_insitu_pot_2d_t 
    650661 
    651662 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traadv.F90

    r14648 r14986  
    1818   USE oce            ! ocean dynamics and active tracers 
    1919   USE dom_oce        ! ocean space and time domain 
    20    ! TEMP: [tiling] This change not necessary after extended haloes development 
     20   ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    2121   USE domtile 
    2222   USE domvvl         ! variable vertical scale factors 
    … …  
    2525   USE traadv_cen     ! centered scheme            (tra_adv_cen  routine) 
    2626   USE traadv_fct     ! FCT      scheme            (tra_adv_fct  routine) 
    27    USE traadv_fct_lf  ! FCT      scheme            (tra_adv_fct  routine - loop fusion version) 
    2827   USE traadv_mus     ! MUSCL    scheme            (tra_adv_mus  routine) 
    29    USE traadv_mus_lf  ! MUSCL    scheme            (tra_adv_mus  routine - loop fusion version) 
    3028   USE traadv_ubs     ! UBS      scheme            (tra_adv_ubs  routine) 
    3129   USE traadv_qck     ! QUICKEST scheme            (tra_adv_qck  routine) 
    … …  
    6159   LOGICAL ::   ln_traadv_qck    ! QUICKEST scheme flag 
    6260 
    63    INTEGER ::   nadv             ! choice of the type of advection scheme 
     61   INTEGER, PUBLIC ::   nadv             ! choice of the type of advection scheme 
    6462   !                             ! associated indices: 
    65    INTEGER, PARAMETER ::   np_NO_adv  = 0   ! no T-S advection 
    66    INTEGER, PARAMETER ::   np_CEN     = 1   ! 2nd/4th order centered scheme 
    67    INTEGER, PARAMETER ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme 
    68    INTEGER, PARAMETER ::   np_MUS     = 3   ! MUSCL scheme 
    69    INTEGER, PARAMETER ::   np_UBS     = 4   ! 3rd order Upstream Biased Scheme 
    70    INTEGER, PARAMETER ::   np_QCK     = 5   ! QUICK scheme 
     63   INTEGER, PARAMETER, PUBLIC ::   np_NO_adv  = 0   ! no T-S advection 
     64   INTEGER, PARAMETER, PUBLIC ::   np_CEN     = 1   ! 2nd/4th order centered scheme 
     65   INTEGER, PARAMETER, PUBLIC ::   np_FCT     = 2   ! 2nd/4th order Flux Corrected Transport scheme 
     66   INTEGER, PARAMETER, PUBLIC ::   np_MUS     = 3   ! MUSCL scheme 
     67   INTEGER, PARAMETER, PUBLIC ::   np_UBS     = 4   ! 3rd order Upstream Biased Scheme 
     68   INTEGER, PARAMETER, PUBLIC ::   np_QCK     = 5   ! QUICK scheme 
    7169 
    7270   !! * Substitutions 
    … …  
    9492      ! 
    9593      INTEGER ::   ji, jj, jk   ! dummy loop index 
    96       ! TEMP: [tiling] This change not necessary and can be A2D(nn_hls) if using XIOS (subdomain support) 
     94      ! TEMP: [tiling] This change not necessary and can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    9795      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zuu, zvv, zww   ! 3D workspace 
    9896      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds 
    99       ! TEMP: [tiling] This change not necessary after extra haloes development 
     97      ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    10098      LOGICAL :: lskip 
    10199      !!---------------------------------------------------------------------- 
    … …  
    105103      lskip = .FALSE. 
    106104 
    107       ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support) 
    108       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     105      ! TEMP: [tiling] These changes not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
     106      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    109107         ALLOCATE( zuu(jpi,jpj,jpk), zvv(jpi,jpj,jpk), zww(jpi,jpj,jpk) ) 
    110108      ENDIF 
    111109 
    112       ! TEMP: [tiling] These changes not necessary after extra haloes development (lbc_lnk removed from tra_adv_*) and if XIOS has subdomain support (ldf_eiv_dia) 
    113       IF( nadv /= np_CEN .OR. (nadv == np_CEN .AND. nn_cen_h == 4) .OR. ln_ldfeiv_dia )  THEN 
    114          IF( ln_tile ) THEN 
    115             IF( ntile == 1 ) THEN 
    116                CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 
    117             ELSE 
    118                lskip = .TRUE. 
    119             ENDIF 
     110      ! TEMP: [tiling] These changes not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
     111      IF( ln_tile .AND. nadv == np_FCT )  THEN 
     112         IF( ntile == 1 ) THEN 
     113            CALL dom_tile_stop( ldhold=.TRUE. ) 
     114         ELSE 
     115            lskip = .TRUE. 
    120116         ENDIF 
    121117      ENDIF 
    … …  
    123119         !                                         !==  effective transport  ==! 
    124120         IF( ln_wave .AND. ln_sdw )  THEN 
    125             DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     121            DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 
    126122               zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * ( uu(ji,jj,jk,Kmm) + usd(ji,jj,jk) ) 
    127123               zvv(ji,jj,jk) = e1v  (ji,jj) * e3v(ji,jj,jk,Kmm) * ( vv(ji,jj,jk,Kmm) + vsd(ji,jj,jk) ) 
    … …  
    129125            END_3D 
    130126         ELSE 
    131             DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     127            DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 
    132128               zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm)               ! eulerian transport only 
    133129               zvv(ji,jj,jk) = e1v  (ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) 
    … …  
    137133         ! 
    138134         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN                                ! add z-tilde and/or vvl corrections 
    139             DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     135            DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 
    140136               zuu(ji,jj,jk) = zuu(ji,jj,jk) + un_td(ji,jj,jk) 
    141137               zvv(ji,jj,jk) = zvv(ji,jj,jk) + vn_td(ji,jj,jk) 
    … …  
    143139         ENDIF 
    144140         ! 
    145          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     141         DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    146142            zuu(ji,jj,jpk) = 0._wp                                                      ! no transport trough the bottom 
    147143            zvv(ji,jj,jpk) = 0._wp 
    … …  
    149145         END_2D 
    150146         ! 
    151          ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support) 
    152147         IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   & 
    153             &              CALL ldf_eiv_trp( kt, nit000, zuu(A2D(nn_hls),:), zvv(A2D(nn_hls),:), zww(A2D(nn_hls),:), & 
    154             &                                'TRA', Kmm, Krhs )   ! add the eiv transport (if necessary) 
    155          ! 
    156          IF( ln_mle    )   CALL tra_mle_trp( kt, nit000, zuu(A2D(nn_hls),:), zvv(A2D(nn_hls),:), zww(A2D(nn_hls),:), & 
    157             &                                'TRA', Kmm       )   ! add the mle transport (if necessary) 
    158          ! 
    159          ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support) 
    160          IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile 
     148            &              CALL ldf_eiv_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm, Krhs )   ! add the eiv transport (if necessary) 
     149         ! 
     150         IF( ln_mle    )   CALL tra_mle_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm       )   ! add the mle transport (if necessary) 
     151         ! 
     152         ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
     153         IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile 
    161154            CALL iom_put( "uocetr_eff", zuu )                                        ! output effective transport 
    162155            CALL iom_put( "vocetr_eff", zvv ) 
    … …  
    164157         ENDIF 
    165158         ! 
    166    !!gm ??? 
    167          ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support) 
     159!!gm ??? 
     160         ! TEMP: [tiling] This copy-in not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    168161         CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) )                                    ! diagnose the effective MSF 
    169    !!gm ??? 
     162!!gm ??? 
    170163         ! 
    171164 
    … …  
    179172         ! 
    180173         CASE ( np_CEN )                                 ! Centered scheme : 2nd / 4th order 
    181             IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kmm), 'T', 1._wp ) 
    182174            CALL tra_adv_cen    ( kt, nit000, 'TRA',         zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v ) 
    183175         CASE ( np_FCT )                                 ! FCT scheme      : 2nd / 4th order 
    184             IF (nn_hls.EQ.2) THEN 
    185                CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1._wp, pts(:,:,:,:,Kmm), 'T', 1._wp) 
    186                CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1._wp, zvv(:,:,:), 'V', -1._wp, zww(:,:,:), 'W', 1._wp) 
    187 #if defined key_loop_fusion 
    188                CALL tra_adv_fct_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 
    189 #else 
    190176               CALL tra_adv_fct    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 
    191 #endif 
    192             ELSE 
    193                CALL tra_adv_fct    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 
    194             END IF 
    195177         CASE ( np_MUS )                                 ! MUSCL 
    196             IF (nn_hls.EQ.2) THEN 
    197                 CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1._wp) 
    198 #if defined key_loop_fusion 
    199                 CALL tra_adv_mus_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 
    200 #else 
    201178                CALL tra_adv_mus    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 
    202 #endif 
    203             ELSE 
    204                 CALL tra_adv_mus    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 
    205             END IF 
    206179         CASE ( np_UBS )                                 ! UBS 
    207             IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1._wp) 
    208180            CALL tra_adv_ubs    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v   ) 
    209181         CASE ( np_QCK )                                 ! QUICKEST 
    210             IF (nn_hls.EQ.2) THEN 
    211                CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1._wp, zvv(:,:,:), 'V', -1._wp) 
    212                CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1._wp) 
    213             END IF 
    214182            CALL tra_adv_qck    ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs ) 
    215183         ! 
    … …  
    226194         ENDIF 
    227195 
    228          ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_adv_*) and if XIOS has subdomain support (ldf_eiv_dia) 
    229          IF( ln_tile .AND. ntile == 0 ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) 
    230  
     196         ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
     197         IF( ln_tile .AND. .NOT. l_istiled ) CALL dom_tile_start( ldhold=.TRUE. ) 
    231198      ENDIF 
    232199      !                                              ! print mean trends (used for debugging) 
    … …  
    234201         &                                  tab3d_2=CASTWP(pts(:,:,:,jp_sal,Krhs)), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    235202 
    236       ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support) 
    237       IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain 
     203      ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
     204      IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only for the full domain 
    238205         DEALLOCATE( zuu, zvv, zww ) 
    239206      ENDIF 
    … …  
    307274        CALL ctl_stop( 'tra_adv_init: FCT scheme, choose 2nd or 4th order' ) 
    308275      ENDIF 
     276      ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
     277      IF( ln_traadv_fct .AND. ln_tile ) THEN 
     278         CALL ctl_warn( 'tra_adv_init: FCT scheme does not yet work with tiling' ) 
     279      ENDIF 
    309280      IF( ln_traadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 )   ) THEN     ! UBS 
    310281        CALL ctl_stop( 'tra_adv_init: UBS scheme, choose 2nd or 4th order' ) 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traadv_cen.F90

    r14644 r14986  
    2323   USE trc_oce        ! share passive tracers/Ocean variables 
    2424   USE lib_mpp        ! MPP library 
     25#if defined key_loop_fusion 
     26   USE traadv_cen_lf  ! centered scheme            (tra_adv_cen  routine - loop fusion version) 
     27#endif 
    2528 
    2629   IMPLICIT NONE 
    … …  
    7275      INTEGER                                  , INTENT(in   ) ::   kn_cen_h        ! =2/4 (2nd or 4th order scheme) 
    7376      INTEGER                                  , INTENT(in   ) ::   kn_cen_v        ! =2/4 (2nd or 4th order scheme) 
    74       ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 
     77      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    7578      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components 
    7679      REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    … …  
    8386      !!---------------------------------------------------------------------- 
    8487      ! 
    85       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     88#if defined key_loop_fusion 
     89      CALL tra_adv_cen_lf    ( kt, nit000, cdtype, pU, pV, pW, Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 
     90#else 
     91      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    8692         IF( kt == kit000 )  THEN 
    8793            IF(lwp) WRITE(numout,*) 
    … …  
    120126               ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
    121127            END_3D 
    122             IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond. 
     128            IF (nn_hls==1) CALL lbc_lnk( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp, ld4only= .TRUE. )   ! Lateral boundary cond. 
    123129            ! 
    124130            DO_3D( nn_hls-1, 0, nn_hls-1, 0, 1, jpkm1 )           ! Horizontal advective fluxes 
    … …  
    132138               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v 
    133139            END_3D 
    134             IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_cen', zwx, 'U', -1._wp , zwy, 'V', -1._wp ) 
     140            IF (nn_hls==1) CALL lbc_lnk( 'traadv_cen', zwx, 'U', -1._wp , zwy, 'V', -1._wp ) 
    135141            ! 
    136142         CASE DEFAULT 
    … …  
    185191      END DO 
    186192      ! 
     193#endif 
    187194   END SUBROUTINE tra_adv_cen 
    188195 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traadv_fct.F90

    r14649 r14986  
    3434   PUBLIC   tra_adv_fct        ! called by traadv.F90 
    3535   PUBLIC   interp_4th_cpt     ! called by traadv_cen.F90 
    36    PUBLIC   tridia_solver      ! called by traadv_fct_lf.F90 
    37    PUBLIC   nonosc             ! called by traadv_fct_lf.F90 - key_agrif 
    3836 
    3937   LOGICAL  ::   l_trd   ! flag to compute trends 
    … …  
    8280      INTEGER                                  , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4) 
    8381      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
    84       ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 
     82      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case 
    8583      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components 
    8684      REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    … …  
    9795      !!---------------------------------------------------------------------- 
    9896      ! 
    99       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     97#if defined key_loop_fusion 
     98      CALL tra_adv_fct_lf ( kt, nit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) 
     99#else 
     100      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    100101         IF( kt == kit000 )  THEN 
    101102            IF(lwp) WRITE(numout,*) 
    … …  
    138139      ! If adaptive vertical advection, check if it is needed on this PE at this time 
    139140      IF( ln_zad_Aimp ) THEN 
    140          IF( MAXVAL( ABS( wi(A2D(nn_hls),:) ) ) > 0._wp ) ll_zAimp = .TRUE. 
     141         IF( MAXVAL( ABS( wi(A2D(1),:) ) ) > 0._wp ) ll_zAimp = .TRUE. 
    141142      END IF 
    142143      ! If active adaptive vertical advection, build tridiagonal matrix 
    … …  
    164165            zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji  ,jj+1,jk,jn,Kbb) ) 
    165166         END_3D 
     167         !                               !* upstream tracer flux in the k direction *! 
     168         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )      ! Interior value ( multiplied by wmask) 
     169            zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 
     170            zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 
     171            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) 
     172         END_3D 
     173         IF( ln_linssh ) THEN               ! top ocean value (only in linear free surface as zwz has been w-masked) 
     174            IF( ln_isfcav ) THEN                        ! top of the ice-shelf cavities and at the ocean surface 
     175               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     176                  zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
     177               END_2D 
     178            ELSE                                        ! no cavities: only at the ocean surface 
     179               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     180                  zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) 
     181               END_2D 
     182            ENDIF 
     183         ENDIF 
     184         ! 
     185         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   !* trend and after field with monotonic scheme 
     186            !                               ! total intermediate advective trends 
     187            ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     188               &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     189               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     190            !                               ! update and guess with monotonic sheme 
     191            pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   & 
     192               &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 
     193            zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 
     194               &                                  / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
     195         END_3D 
     196 
     197         IF ( ll_zAimp ) THEN 
     198            CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) 
     199            ! 
     200            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 
     201            DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Interior value ( multiplied by wmask) 
     202               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     203               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     204               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) 
     205               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 
     206            END_3D 
     207            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     208               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
     209                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     210            END_3D 
     211            ! 
     212         END IF 
     213         ! 
     214         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
     215            ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:) 
     216         END IF 
     217         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
     218         IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:) 
     219         ! 
     220         !        !==  anti-diffusive flux : high order minus low order  ==! 
     221         ! 
     222         SELECT CASE( kn_fct_h )    !* horizontal anti-diffusive fluxes 
     223         ! 
     224         CASE(  2  )                   !- 2nd order centered 
     225            DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 
     226               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) 
     227               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) 
     228            END_3D 
     229            ! 
     230         CASE(  4  )                   !- 4th order centered 
     231            zltu(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero 
     232            zltv(:,:,jpk) = 0._wp 
     233            DO jk = 1, jpkm1                 ! Laplacian 
     234               DO_2D( 1, 0, 1, 0 )                 ! 1st derivative (gradient) 
     235                  ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 
     236                  ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
     237               END_2D 
     238               DO_2D( 0, 0, 0, 0 )                 ! 2nd derivative * 1/ 6 
     239                  zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6 
     240                  zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6 
     241               END_2D 
     242            END DO 
     243            ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility 
     244            CALL lbc_lnk( 'traadv_fct', zltu, 'T', -1.0_wp , zltv, 'T', -1.0_wp, ld4only= .TRUE. )   ! Lateral boundary cond. (unchanged sgn) 
     245            ! 
     246            DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 
     247               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 
     248               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) 
     249               !                                                        ! C4 minus upstream advective fluxes 
     250               ! round brackets added to fix the order of floating point operations 
     251               ! needed to ensure halo 1 - halo 2 compatibility 
     252               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu(ji,jj,jk) - zltu(ji+1,jj,jk)   & 
     253                             &                                     )                                     & ! bracket for halo 1 - halo 2 compatibility 
     254                             &                          ) - zwx(ji,jj,jk) 
     255               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv(ji,jj,jk) - zltv(ji,jj+1,jk)   & 
     256                             &                                     )                                     & ! bracket for halo 1 - halo 2 compatibility 
     257                             &                          ) - zwy(ji,jj,jk) 
     258            END_3D 
     259            ! 
     260         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested 
     261            ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero 
     262            ztv(:,:,jpk) = 0._wp 
     263            DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )    ! 1st derivative (gradient) 
     264               ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 
     265               ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
     266            END_3D 
     267            IF (nn_hls==1) CALL lbc_lnk( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp, ld4only= .TRUE. )   ! Lateral boundary cond. (unchanged sgn) 
     268            ! 
     269            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Horizontal advective fluxes 
     270               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) 
     271               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) 
     272               !                                                  ! C4 interpolation of T at u- & v-points (x2) 
     273               zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) ) 
     274               zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) ) 
     275               !                                                  ! C4 minus upstream advective fluxes 
     276               zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 
     277               zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 
     278            END_3D 
     279            IF (nn_hls==2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
     280            ! 
     281         END SELECT 
     282         ! 
     283         SELECT CASE( kn_fct_v )    !* vertical anti-diffusive fluxes (w-masked interior values) 
     284         ! 
     285         CASE(  2  )                   !- 2nd order centered 
     286            DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
     287               zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   & 
     288                  &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk) 
     289            END_3D 
     290            ! 
     291         CASE(  4  )                   !- 4th order COMPACT 
     292            CALL interp_4th_cpt( CASTWP(pt(:,:,:,jn,Kmm)) , ztw )   ! zwt = COMPACT interpolation of T at w-point 
     293            DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
     294               zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
     295            END_3D 
     296            ! 
     297         END SELECT 
     298         IF( ln_linssh ) THEN    ! top ocean value: high order = upstream  ==>>  zwz=0 
     299            zwz(:,:,1) = 0._wp   ! only ocean surface as interior zwz values have been w-masked 
     300         ENDIF 
     301         ! 
     302         IF (nn_hls==1) THEN 
     303            CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp ) 
     304            CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp ) 
     305         ELSE 
     306            CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 
     307         END IF 
     308         ! 
     309         IF ( ll_zAimp ) THEN 
     310            DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )    !* trend and after field with monotonic scheme 
     311               !                                                ! total intermediate advective trends 
     312               ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     313                  &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     314                  &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     315               ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
     316            END_3D 
     317            ! 
     318            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 
     319            ! 
     320            DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Interior value ( multiplied by wmask) 
     321               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
     322               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
     323               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) 
     324            END_3D 
     325         END IF 
     326         ! 
     327         !        !==  monotonicity algorithm  ==! 
     328         ! 
     329         CALL nonosc( Kmm, CASTWP(pt(:,:,:,jn,Kbb)), zwx, zwy, zwz, zwi, p2dt ) 
     330         ! 
     331         !        !==  final trend with corrected fluxes  ==! 
     332         ! 
     333         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     334            ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
     335               &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     336               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
     337            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) 
     338            zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
     339         END_3D 
     340         ! 
     341         IF ( ll_zAimp ) THEN 
     342            ! 
     343            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 
     344            DO_3D( 0, 0, 0, 0, 2, jpkm1 )      ! Interior value ( multiplied by wmask) 
     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               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) 
     348               zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 
     349            END_3D 
     350            DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
     351               pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji  ,jj  ,jk+1) ) & 
     352                  &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     353            END_3D 
     354         END IF 
     355         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport 
     356            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< add anti-diffusive fluxes 
     357            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  !     to upstream fluxes 
     358            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! 
     359            ! 
     360            IF( l_trd ) THEN              ! trend diagnostics 
     361               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, CASTWP(pt(:,:,:,jn,Kmm)) ) 
     362               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, CASTWP(pt(:,:,:,jn,Kmm)) ) 
     363               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, CASTWP(pt(:,:,:,jn,Kmm)) ) 
     364            ENDIF 
     365            !                             ! heat/salt transport 
     366            IF( l_hst )   CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 
     367            ! 
     368         ENDIF 
     369         IF( l_ptr ) THEN              ! "Poleward" transports 
     370            zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< add anti-diffusive fluxes 
     371            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 
     372         ENDIF 
     373         ! 
     374      END DO                     ! end of tracer loop 
     375      ! 
     376      IF ( ll_zAimp ) THEN 
     377         DEALLOCATE( zwdia, zwinf, zwsup ) 
     378      ENDIF 
     379      IF( l_trd .OR. l_hst ) THEN 
     380         DEALLOCATE( ztrdx, ztrdy, ztrdz ) 
     381      ENDIF 
     382      IF( l_ptr ) THEN 
     383         DEALLOCATE( zptry ) 
     384      ENDIF 
     385      ! 
     386#endif 
     387   END SUBROUTINE tra_adv_fct 
     388 
     389 
     390   SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt ) 
     391      !!--------------------------------------------------------------------- 
     392      !!                    ***  ROUTINE nonosc  *** 
     393      !! 
     394      !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
     395      !!       scheme and the before field by a nonoscillatory algorithm 
     396      !! 
     397      !! **  Method  :   ... ??? 
     398      !!       warning : pbef and paft must be masked, but the boundaries 
     399      !!       conditions on the fluxes are not necessary zalezak (1979) 
     400      !!       drange (1995) multi-dimensional forward-in-time and upstream- 
     401      !!       in-space based differencing for fluid 
     402      !!---------------------------------------------------------------------- 
     403      INTEGER                         , INTENT(in   ) ::   Kmm             ! time level index 
     404      REAL(wp)                        , INTENT(in   ) ::   p2dt            ! tracer time-step 
     405      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pbef            ! before field 
     406      REAL(wp), DIMENSION(A2D(nn_hls)    ,jpk), INTENT(in   ) ::   paft            ! after field 
     407      REAL(dp), DIMENSION(A2D(nn_hls)    ,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions 
     408      ! 
     409      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     410      INTEGER  ::   ikm1         ! local integer 
     411      REAL(dp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars 
     412      REAL(dp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
     413      REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo 
     414      !!---------------------------------------------------------------------- 
     415      ! 
     416      zbig  = 1.e+40_dp 
     417      zrtrn = 1.e-15_dp 
     418      zbetup(:,:,:) = 0._dp   ;   zbetdo(:,:,:) = 0._dp 
     419 
     420      ! Search local extrema 
     421      ! -------------------- 
     422      ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
     423      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 
     424         zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ),   & 
     425            &                  paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) )  ) 
     426         zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ),   & 
     427            &                  paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) )  ) 
     428      END_3D 
     429 
     430      DO jk = 1, jpkm1 
     431         ikm1 = MAX(jk-1,1) 
     432         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     433 
     434            ! search maximum in neighbourhood 
     435            zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 
     436               &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   & 
     437               &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   & 
     438               &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  ) 
     439 
     440            ! search minimum in neighbourhood 
     441            zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   & 
     442               &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   & 
     443               &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   & 
     444               &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  ) 
     445 
     446            ! positive part of the flux 
     447            zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   & 
     448               & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   & 
     449               & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) ) 
     450 
     451            ! negative part of the flux 
     452            zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   & 
     453               & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   & 
     454               & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
     455 
     456            ! up & down beta terms 
     457            zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt 
     458            zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 
     459            zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt 
     460         END_2D 
     461      END DO 
     462      IF (nn_hls==1) CALL lbc_lnk( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp, ld4only= .TRUE. )   ! lateral boundary cond. (unchanged sign) 
     463 
     464      ! 3. monotonic flux in the i & j direction (paa & pbb) 
     465      ! ---------------------------------------- 
     466      DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
     467         zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
     468         zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
     469         zcu =       ( 0.5  + SIGN( 0.5_wp , CASTWP(paa(ji,jj,jk)) ) ) 
     470         paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
     471 
     472         zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
     473         zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
     474         zcv =       ( 0.5  + SIGN( 0.5_wp , CASTWP(pbb(ji,jj,jk)) ) ) 
     475         pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
     476 
     477      ! monotonic flux in the k direction, i.e. pcc 
     478      ! ------------------------------------------- 
     479         za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 
     480         zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 
     481         zc =       ( 0.5  + SIGN( 0.5_wp , CASTWP(pcc(ji,jj,jk+1)) ) ) 
     482         pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 
     483      END_3D 
     484      ! 
     485   END SUBROUTINE nonosc 
     486 
     487 
     488   SUBROUTINE interp_4th_cpt_org( pt_in, pt_out ) 
     489      !!---------------------------------------------------------------------- 
     490      !!                  ***  ROUTINE interp_4th_cpt_org  *** 
     491      !! 
     492      !! **  Purpose :   Compute the interpolation of tracer at w-point 
     493      !! 
     494      !! **  Method  :   4th order compact interpolation 
     495      !!---------------------------------------------------------------------- 
     496      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! now tracer fields 
     497      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! now tracer field interpolated at w-pts 
     498      ! 
     499      INTEGER :: ji, jj, jk   ! dummy loop integers 
     500      REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 
     501      !!---------------------------------------------------------------------- 
     502 
     503      DO_3D( 1, 1, 1, 1, 3, jpkm1 )       !==  build the three diagonal matrix  ==! 
     504         zwd (ji,jj,jk) = 4._wp 
     505         zwi (ji,jj,jk) = 1._wp 
     506         zws (ji,jj,jk) = 1._wp 
     507         zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
     508         ! 
     509         IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom 
     510            zwd (ji,jj,jk) = 1._wp 
     511            zwi (ji,jj,jk) = 0._wp 
     512            zws (ji,jj,jk) = 0._wp 
     513            zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
     514         ENDIF 
     515      END_3D 
     516      ! 
     517      jk = 2                                    ! Switch to second order centered at top 
     518      DO_2D( 1, 1, 1, 1 ) 
     519         zwd (ji,jj,jk) = 1._wp 
     520         zwi (ji,jj,jk) = 0._wp 
     521         zws (ji,jj,jk) = 0._wp 
     522         zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
     523      END_2D 
     524      ! 
     525      !                       !==  tridiagonal solve  ==! 
     526      DO_2D( 1, 1, 1, 1 )           ! first recurrence 
     527         zwt(ji,jj,2) = zwd(ji,jj,2) 
     528      END_2D 
     529      DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 
     530         zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
     531      END_3D 
     532      ! 
     533      DO_2D( 1, 1, 1, 1 )           ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
     534         pt_out(ji,jj,2) = zwrm(ji,jj,2) 
     535      END_2D 
     536      DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 
     537         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 
     538      END_3D 
     539 
     540      DO_2D( 1, 1, 1, 1 )           ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 
     541         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
     542      END_2D 
     543      DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 ) 
     544         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
     545      END_3D 
     546      ! 
     547   END SUBROUTINE interp_4th_cpt_org 
     548 
     549 
     550   SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 
     551      !!---------------------------------------------------------------------- 
     552      !!                  ***  ROUTINE interp_4th_cpt  *** 
     553      !! 
     554      !! **  Purpose :   Compute the interpolation of tracer at w-point 
     555      !! 
     556      !! **  Method  :   4th order compact interpolation 
     557      !!---------------------------------------------------------------------- 
     558      REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! field at t-point 
     559      REAL(wp),DIMENSION(A2D(nn_hls)    ,jpk), INTENT(  out) ::   pt_out   ! field interpolated at w-point 
     560      ! 
     561      INTEGER ::   ji, jj, jk   ! dummy loop integers 
     562      INTEGER ::   ikt, ikb     ! local integers 
     563      REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwd, zwi, zws, zwrm, zwt 
     564      !!---------------------------------------------------------------------- 
     565      ! 
     566      !                      !==  build the three diagonal matrix & the RHS  ==! 
     567      ! 
     568      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 )    ! interior (from jk=3 to jpk-1) 
     569         zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal 
     570         zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal 
     571         zws (ji,jj,jk) =         wmask(ji,jj,jk)                         ! upper diagonal 
     572         zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk)                     &   ! RHS 
     573            &           *       ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 
     574      END_3D 
     575      ! 
     576!!gm 
     577!      SELECT CASE( kbc )               !* boundary condition 
     578!      CASE( np_NH   )   ! Neumann homogeneous at top & bottom 
     579!      CASE( np_CEN2 )   ! 2nd order centered  at top & bottom 
     580!      END SELECT 
     581!!gm 
     582      ! 
     583      IF ( ln_isfcav ) THEN            ! set level two values which may not be set in ISF case 
     584         zwd(:,:,2) = 1._wp  ;  zwi(:,:,2) = 0._wp  ;  zws(:,:,2) = 0._wp  ;  zwrm(:,:,2) = 0._wp 
     585      END IF 
     586      ! 
     587      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )              ! 2nd order centered at top & bottom 
     588         ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point 
     589         ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point 
     590         ! 
     591         zwd (ji,jj,ikt) = 1._wp          ! top 
     592         zwi (ji,jj,ikt) = 0._wp 
     593         zws (ji,jj,ikt) = 0._wp 
     594         zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) ) 
     595         ! 
     596         zwd (ji,jj,ikb) = 1._wp          ! bottom 
     597         zwi (ji,jj,ikb) = 0._wp 
     598         zws (ji,jj,ikb) = 0._wp 
     599         zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 
     600      END_2D 
     601      ! 
     602      !                       !==  tridiagonal solver  ==! 
     603      ! 
     604      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )           !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1 
     605         zwt(ji,jj,2) = zwd(ji,jj,2) 
     606      END_2D 
     607      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 
     608         zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
     609      END_3D 
     610      ! 
     611      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
     612         pt_out(ji,jj,2) = zwrm(ji,jj,2) 
     613      END_2D 
     614      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 
     615         pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 
     616      END_3D 
     617 
     618      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk 
     619         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
     620      END_2D 
     621      DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 
     622         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
     623      END_3D 
     624      ! 
     625   END SUBROUTINE interp_4th_cpt 
     626 
     627 
     628   SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev ) 
     629      !!---------------------------------------------------------------------- 
     630      !!                  ***  ROUTINE tridia_solver  *** 
     631      !! 
     632      !! **  Purpose :   solve a symmetric 3diagonal system 
     633      !! 
     634      !! **  Method  :   solve M.t_out = RHS(t)  where M is a tri diagonal matrix ( jpk*jpk ) 
     635      !! 
     636      !!             ( D_1 U_1  0   0   0  )( t_1 )   ( RHS_1 ) 
     637      !!             ( L_2 D_2 U_2  0   0  )( t_2 )   ( RHS_2 ) 
     638      !!             (  0  L_3 D_3 U_3  0  )( t_3 ) = ( RHS_3 ) 
     639      !!             (        ...          )( ... )   ( ...  ) 
     640      !!             (  0   0   0  L_k D_k )( t_k )   ( RHS_k ) 
     641      !! 
     642      !!        M is decomposed in the product of an upper and lower triangular matrix. 
     643      !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL 
     644      !!        (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 
     645      !!        The solution is pta. 
     646      !!        The 3d array zwt is used as a work space array. 
     647      !!---------------------------------------------------------------------- 
     648      REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pD, pU, pL    ! 3-diagonal matrix 
     649      REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pRHS          ! Right-Hand-Side 
     650      REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(  out) ::   pt_out        !!gm field at level=F(klev) 
     651      INTEGER                    , INTENT(in   ) ::   klev          ! =1 pt_out at w-level 
     652      !                                                             ! =0 pt at t-level 
     653      INTEGER ::   ji, jj, jk   ! dummy loop integers 
     654      INTEGER ::   kstart       ! local indices 
     655      REAL(wp),DIMENSION(A2D(nn_hls),jpk) ::   zwt   ! 3D work array 
     656      !!---------------------------------------------------------------------- 
     657      ! 
     658      kstart =  1  + klev 
     659      ! 
     660      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                         !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1 
     661         zwt(ji,jj,kstart) = pD(ji,jj,kstart) 
     662      END_2D 
     663      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 
     664         zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
     665      END_3D 
     666      ! 
     667      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                        !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
     668         pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 
     669      END_2D 
     670      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 
     671         pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 
     672      END_3D 
     673 
     674      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                       !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk 
     675         pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
     676      END_2D 
     677      DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, kstart, -1 ) 
     678         pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
     679      END_3D 
     680      ! 
     681   END SUBROUTINE tridia_solver 
     682 
     683#if defined key_loop_fusion 
     684#define tracer_flux_i(out,zfp,zfm,ji,jj,jk) \ 
     685        zfp = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) ; \ 
     686        zfm = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) ) ; \ 
     687        out = 0.5 * ( zfp * pt(ji,jj,jk,jn,Kbb) + zfm * pt(ji+1,jj,jk,jn,Kbb) ) 
     688 
     689#define tracer_flux_j(out,zfp,zfm,ji,jj,jk) \ 
     690        zfp = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) ) ; \ 
     691        zfm = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) ) ; \ 
     692        out = 0.5 * ( zfp * pt(ji,jj,jk,jn,Kbb) + zfm * pt(ji,jj+1,jk,jn,Kbb) ) 
     693 
     694   SUBROUTINE tra_adv_fct_lf( kt, kit000, cdtype, p2dt, pU, pV, pW,       & 
     695      &                    Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) 
     696      !!---------------------------------------------------------------------- 
     697      !!                  ***  ROUTINE tra_adv_fct  *** 
     698      !! 
     699      !! **  Purpose :   Compute the now trend due to total advection of tracers 
     700      !!               and add it to the general trend of tracer equations 
     701      !! 
     702      !! **  Method  : - 2nd or 4th FCT scheme on the horizontal direction 
     703      !!               (choice through the value of kn_fct) 
     704      !!               - on the vertical the 4th order is a compact scheme 
     705      !!               - corrected flux (monotonic correction) 
     706      !! 
     707      !! ** Action : - update pt(:,:,:,:,Krhs)  with the now advective tracer trends 
     708      !!             - send trends to trdtra module for further diagnostics (l_trdtra=T) 
     709      !!             - poleward advective heat and salt transport (ln_diaptr=T) 
     710      !!---------------------------------------------------------------------- 
     711      INTEGER                                  , INTENT(in   ) ::   kt              ! ocean time-step index 
     712      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs  ! ocean time level indices 
     713      INTEGER                                  , INTENT(in   ) ::   kit000          ! first time step index 
     714      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator) 
     715      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers 
     716      INTEGER                                  , INTENT(in   ) ::   kn_fct_h        ! order of the FCT scheme (=2 or 4) 
     717      INTEGER                                  , INTENT(in   ) ::   kn_fct_v        ! order of the FCT scheme (=2 or 4) 
     718      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
     719      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components 
     720      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
     721      ! 
     722      INTEGER  ::   ji, jj, jk, jn                           ! dummy loop indices 
     723      REAL(wp) ::   ztra                                     ! local scalar 
     724      REAL(wp) ::   zwx_im1, zfp_ui, zfp_ui_m1, zfp_vj, zfp_vj_m1, zfp_wk, zC2t_u, zC4t_u   !   -      - 
     725      REAL(wp) ::   zwy_jm1, zfm_ui, zfm_ui_m1, zfm_vj, zfm_vj_m1, zfm_wk, zC2t_v, zC4t_v   !   -      - 
     726      REAL(wp) ::   ztu, ztv, ztu_im1, ztu_ip1, ztv_jm1, ztv_jp1 
     727      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx_3d, zwy_3d, zwz, ztw, zltu_3d, zltv_3d 
     728      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdx, ztrdy, ztrdz, zptry 
     729      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zwinf, zwdia, zwsup 
     730      LOGICAL  ::   ll_zAimp                                 ! flag to apply adaptive implicit vertical advection 
     731      !!---------------------------------------------------------------------- 
     732      ! 
     733      IF( kt == kit000 )  THEN 
     734         IF(lwp) WRITE(numout,*) 
     735         IF(lwp) WRITE(numout,*) 'tra_adv_fct_lf : FCT advection scheme on ', cdtype 
     736         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
     737      ENDIF 
     738      !! -- init to 0 
     739      zwx_3d(:,:,:) = 0._wp 
     740      zwy_3d(:,:,:) = 0._wp 
     741      zwz(:,:,:) = 0._wp 
     742      zwi(:,:,:) = 0._wp 
     743      ! 
     744      l_trd = .FALSE.            ! set local switches 
     745      l_hst = .FALSE. 
     746      l_ptr = .FALSE. 
     747      ll_zAimp = .FALSE. 
     748      IF( ( cdtype == 'TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )      l_trd = .TRUE. 
     749      IF(   cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) )    l_ptr = .TRUE. 
     750      IF(   cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR.  & 
     751         &                         iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE. 
     752      ! 
     753      IF( l_trd .OR. l_hst )  THEN 
     754         ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) ) 
     755         ztrdx(:,:,:) = 0._wp   ;    ztrdy(:,:,:) = 0._wp   ;   ztrdz(:,:,:) = 0._wp 
     756      ENDIF 
     757      ! 
     758      IF( l_ptr ) THEN 
     759         ALLOCATE( zptry(jpi,jpj,jpk) ) 
     760         zptry(:,:,:) = 0._wp 
     761      ENDIF 
     762      ! 
     763      ! If adaptive vertical advection, check if it is needed on this PE at this time 
     764      IF( ln_zad_Aimp ) THEN 
     765         IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE. 
     766      END IF 
     767      ! If active adaptive vertical advection, build tridiagonal matrix 
     768      IF( ll_zAimp ) THEN 
     769         ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 
     770         DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 
     771            zwdia(ji,jj,jk) =  1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )   & 
     772            &                               / e3t(ji,jj,jk,Krhs) 
     773            zwinf(ji,jj,jk) =  p2dt * MIN( wi(ji,jj,jk  ) , 0._wp ) / e3t(ji,jj,jk,Krhs) 
     774            zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) 
     775         END_3D 
     776      END IF 
     777      ! 
     778      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
     779         ! 
     780         !        !==  upstream advection with initial mass fluxes & intermediate update  ==! 
    166781         !                               !* upstream tracer flux in the k direction *! 
    167782         DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! Interior value ( multiplied by wmask) 
    … …  
    182797         ENDIF 
    183798         ! 
    184          DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   !* trend and after field with monotonic scheme 
    185             !                               ! total intermediate advective trends 
    186             ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    187                &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
    188                &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
    189             !                               ! update and guess with monotonic sheme 
    190             pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   & 
    191                &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 
    192             zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 
    193                &                                  / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
    194          END_3D 
     799         !                    !* upstream tracer flux in the i and j direction 
     800         DO jk = 1, jpkm1 
     801            DO jj = 1, jpj-1 
     802               tracer_flux_i(zwx_3d(1,jj,jk),zfp_ui,zfm_ui,1,jj,jk) 
     803               tracer_flux_j(zwy_3d(1,jj,jk),zfp_vj,zfm_vj,1,jj,jk) 
     804            END DO 
     805            DO ji = 1, jpi-1 
     806               tracer_flux_i(zwx_3d(ji,1,jk),zfp_ui,zfm_ui,ji,1,jk) 
     807               tracer_flux_j(zwy_3d(ji,1,jk),zfp_vj,zfm_vj,ji,1,jk) 
     808            END DO 
     809            DO_2D( 1, 1, 1, 1 ) 
     810               tracer_flux_i(zwx_3d(ji,jj,jk),zfp_ui,zfm_ui,ji,jj,jk) 
     811               tracer_flux_i(zwx_im1,zfp_ui_m1,zfm_ui_m1,ji-1,jj,jk) 
     812               tracer_flux_j(zwy_3d(ji,jj,jk),zfp_vj,zfm_vj,ji,jj,jk) 
     813               tracer_flux_j(zwy_jm1,zfp_vj_m1,zfm_vj_m1,ji,jj-1,jk) 
     814               ztra = - ( zwx_3d(ji,jj,jk) - zwx_im1 + zwy_3d(ji,jj,jk) - zwy_jm1 + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) 
     815               !                               ! update and guess with monotonic sheme 
     816               pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   & 
     817                  &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 
     818               zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 
     819                  &                                  / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
     820            END_2D 
     821         END DO 
    195822 
    196823         IF ( ll_zAimp ) THEN 
    … …  
    198825            ! 
    199826            ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 
    200             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Interior value ( multiplied by wmask) 
     827            DO_3D( 1, 1, 1, 1, 2, jpkm1 )       ! Interior value ( multiplied by wmask) 
    201828               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
    202829               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
    … …  
    212839         ! 
    213840         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
    214             ztrdx(:,:,:) = zwx(:,:,:)   ;   ztrdy(:,:,:) = zwy(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:) 
     841            ztrdx(:,:,:) = zwx_3d(:,:,:)   ;   ztrdy(:,:,:) = zwy_3d(:,:,:)   ;   ztrdz(:,:,:) = zwz(:,:,:) 
    215842         END IF 
    216843         !                             ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    217          IF( l_ptr )   zptry(:,:,:) = zwy(:,:,:) 
     844         IF( l_ptr )   zptry(:,:,:) = zwy_3d(:,:,:) 
    218845         ! 
    219846         !        !==  anti-diffusive flux : high order minus low order  ==! 
    … …  
    222849         ! 
    223850         CASE(  2  )                   !- 2nd order centered 
    224             DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 
    225                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) 
    226                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) 
     851            DO_3D( 2, 1, 2, 1, 1, jpkm1 ) 
     852               zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx_3d(ji,jj,jk) 
     853               zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy_3d(ji,jj,jk) 
    227854            END_3D 
    228855            ! 
    229856         CASE(  4  )                   !- 4th order centered 
    230             zltu(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero 
    231             zltv(:,:,jpk) = 0._wp 
    232             DO jk = 1, jpkm1                 ! Laplacian 
    233                DO_2D( 1, 0, 1, 0 )                 ! 1st derivative (gradient) 
    234                   ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 
    235                   ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
    236                END_2D 
    237                DO_2D( 0, 0, 0, 0 )                 ! 2nd derivative * 1/ 6 
    238                   zltu(ji,jj,jk) = (  ztu(ji,jj,jk) + ztu(ji-1,jj,jk)  ) * r1_6 
    239                   zltv(ji,jj,jk) = (  ztv(ji,jj,jk) + ztv(ji,jj-1,jk)  ) * r1_6 
     857            zltu_3d(:,:,jpk) = 0._wp            ! Bottom value : flux set to zero 
     858            zltv_3d(:,:,jpk) = 0._wp 
     859            !                          ! Laplacian 
     860            DO_3D( 0, 0, 0, 0, 1, jpkm1 )                 ! 2nd derivative * 1/ 6 
     861                  !             ! 1st derivative (gradient) 
     862                  ztu = ( pt(ji+1,jj,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 
     863                  ztu_im1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk) 
     864                  ztv = ( pt(ji,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
     865                  ztv_jm1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk) 
     866                  !             ! 2nd derivative * 1/ 6 
     867                  zltu_3d(ji,jj,jk) = (  ztu + ztu_im1  ) * r1_6 
     868                  zltv_3d(ji,jj,jk) = (  ztv + ztv_jm1  ) * r1_6 
    240869               END_2D 
    241870            END DO 
    242             CALL lbc_lnk( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
    243             ! 
    244             DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
     871            ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility 
     872            CALL lbc_lnk( 'traadv_fct', zltu_3d, 'T', -1.0_wp , zltv_3d, 'T', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
     873            ! 
     874            DO_3D( 2, 1, 2, 1, 1, jpkm1 ) 
    245875               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 
    246876               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) 
    247877               !                                                        ! C4 minus upstream advective fluxes 
    248                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) 
    249                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) 
    250             END_3D 
    251             IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp, zwy, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
     878               ! round brackets added to fix the order of floating point operations 
     879               ! needed to ensure halo 1 - halo 2 compatibility 
     880               zwx_3d(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu_3d(ji,jj,jk) - zltu_3d(ji+1,jj,jk)   & 
     881                             &                                        )                                           & ! bracket for halo 1 - halo 2 compatibility 
     882                             &                             ) - zwx_3d(ji,jj,jk) 
     883               zwy_3d(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv_3d(ji,jj,jk) - zltv_3d(ji,jj+1,jk)   & 
     884                             &                                        )                                           & ! bracket for halo 1 - halo 2 compatibility 
     885                             &                             ) - zwy_3d(ji,jj,jk) 
     886            END_3D 
    252887            ! 
    253888         CASE(  41 )                   !- 4th order centered       ==>>   !!gm coding attempt   need to be tested 
    254             ztu(:,:,jpk) = 0._wp             ! Bottom value : flux set to zero 
    255             ztv(:,:,jpk) = 0._wp 
    256             DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )    ! 1st derivative (gradient) 
    257                ztu(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 
    258                ztv(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 
    259             END_3D 
    260             IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
    261             ! 
    262889            DO_3D( 0, 0, 0, 0, 1, jpkm1 )    ! Horizontal advective fluxes 
     890               ztu_im1 = ( pt(ji  ,jj  ,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk) 
     891               ztu_ip1 = ( pt(ji+2,jj  ,jk,jn,Kmm) - pt(ji+1,jj,jk,jn,Kmm) ) * umask(ji+1,jj,jk) 
     892 
     893               ztv_jm1 = ( pt(ji,jj  ,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk) 
     894               ztv_jp1 = ( pt(ji,jj+2,jk,jn,Kmm) - pt(ji,jj+1,jk,jn,Kmm) ) * vmask(ji,jj+1,jk) 
    263895               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) 
    264896               zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji  ,jj+1,jk,jn,Kmm) 
    265897               !                                                  ! C4 interpolation of T at u- & v-points (x2) 
    266                zC4t_u =  zC2t_u + r1_6 * ( ztu(ji-1,jj  ,jk) - ztu(ji+1,jj  ,jk) ) 
    267                zC4t_v =  zC2t_v + r1_6 * ( ztv(ji  ,jj-1,jk) - ztv(ji  ,jj+1,jk) ) 
     898               zC4t_u =  zC2t_u + r1_6 * ( ztu_im1 - ztu_ip1 ) 
     899               zC4t_v =  zC2t_v + r1_6 * ( ztv_jm1 - ztv_jp1 ) 
    268900               !                                                  ! C4 minus upstream advective fluxes 
    269                zwx(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 
    270                zwy(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 
    271             END_3D 
    272             IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
     901               zwx_3d(ji,jj,jk) =  0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx_3d(ji,jj,jk) 
     902               zwy_3d(ji,jj,jk) =  0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy_3d(ji,jj,jk) 
     903            END_3D 
     904            CALL lbc_lnk( 'traadv_fct', zwx_3d, 'U', -1.0_wp , zwy_3d, 'V', -1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
    273905            ! 
    274906         END SELECT 
    … …  
    277909         ! 
    278910         CASE(  2  )                   !- 2nd order centered 
    279             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
     911            DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 
    280912               zwz(ji,jj,jk) =  (  pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   & 
    281913                  &              - zwz(ji,jj,jk)  ) * wmask(ji,jj,jk) 
    … …  
    283915            ! 
    284916         CASE(  4  )                   !- 4th order COMPACT 
    285             CALL interp_4th_cpt( CASTWP(pt(:,:,:,jn,Kmm)) , ztw )   ! zwt = COMPACT interpolation of T at w-point 
    286             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
     917            CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw )   ! zwt = COMPACT interpolation of T at w-point 
     918            DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 
    287919               zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 
    288920            END_3D 
    … …  
    293925         ENDIF 
    294926         ! 
    295          IF (nn_hls.EQ.1) THEN 
    296             CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp ) 
    297             CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp ) 
    298          ELSE 
    299             CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 
    300          END IF 
     927         CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 
    301928         ! 
    302929         IF ( ll_zAimp ) THEN 
    303             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )    !* trend and after field with monotonic scheme 
     930            DO_3D( 1, 1, 1, 1, 1, jpkm1 )    !* trend and after field with monotonic scheme 
    304931               !                                                ! total intermediate advective trends 
    305                ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    306                   &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     932               ztra = - (  zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj  ,jk  )   & 
     933                  &      + zwy_3d(ji,jj,jk) - zwy_3d(ji  ,jj-1,jk  )   & 
    307934                  &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
    308935               ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 
    … …  
    311938            CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 
    312939            ! 
    313             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! Interior value ( multiplied by wmask) 
     940            DO_3D( 1, 1, 1, 1, 2, jpkm1 )       ! Interior value ( multiplied by wmask) 
    314941               zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 
    315942               zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 
    … …  
    320947         !        !==  monotonicity algorithm  ==! 
    321948         ! 
    322          CALL nonosc( Kmm, CASTWP(pt(:,:,:,jn,Kbb)), zwx, zwy, zwz, zwi, p2dt ) 
     949         CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt ) 
    323950         ! 
    324951         !        !==  final trend with corrected fluxes  ==! 
    325952         ! 
    326953         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    327             ztra = - (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   & 
    328                &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   & 
     954            ztra = - (  zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj  ,jk  )   & 
     955               &      + zwy_3d(ji,jj,jk) - zwy_3d(ji  ,jj-1,jk  )   & 
    329956               &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj) 
    330957            pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) 
    … …  
    346973            END_3D 
    347974         END IF 
     975         ! NOT TESTED - NEED l_trd OR l_hst TRUE 
    348976         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport 
    349             ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:)  ! <<< add anti-diffusive fluxes 
    350             ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:)  !     to upstream fluxes 
     977            ztrdx(:,:,:) = ztrdx(:,:,:) + zwx_3d(:,:,:)  ! <<< add anti-diffusive fluxes 
     978            ztrdy(:,:,:) = ztrdy(:,:,:) + zwy_3d(:,:,:)  !     to upstream fluxes 
    351979            ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:)  ! 
    352980            ! 
    353981            IF( l_trd ) THEN              ! trend diagnostics 
    354                CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, CASTWP(pt(:,:,:,jn,Kmm)) ) 
    355                CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, CASTWP(pt(:,:,:,jn,Kmm)) ) 
    356                CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, CASTWP(pt(:,:,:,jn,Kmm)) ) 
     982               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) ) 
     983               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) ) 
     984               CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) ) 
    357985            ENDIF 
    358986            !                             ! heat/salt transport 
    … …  
    360988            ! 
    361989         ENDIF 
     990         ! NOT TESTED - NEED l_ptr TRUE 
    362991         IF( l_ptr ) THEN              ! "Poleward" transports 
    363             zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:)  ! <<< add anti-diffusive fluxes 
     992            zptry(:,:,:) = zptry(:,:,:) + zwy_3d(:,:,:)  ! <<< add anti-diffusive fluxes 
    364993            CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 
    365994         ENDIF 
    … …  
    3771006      ENDIF 
    3781007      ! 
    379    END SUBROUTINE tra_adv_fct 
    380  
    381  
    382    SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt ) 
    383       !!--------------------------------------------------------------------- 
    384       !!                    ***  ROUTINE nonosc  *** 
    385       !! 
    386       !! **  Purpose :   compute monotonic tracer fluxes from the upstream 
    387       !!       scheme and the before field by a nonoscillatory algorithm 
    388       !! 
    389       !! **  Method  :   ... ??? 
    390       !!       warning : pbef and paft must be masked, but the boundaries 
    391       !!       conditions on the fluxes are not necessary zalezak (1979) 
    392       !!       drange (1995) multi-dimensional forward-in-time and upstream- 
    393       !!       in-space based differencing for fluid 
    394       !!---------------------------------------------------------------------- 
    395       INTEGER                         , INTENT(in   ) ::   Kmm             ! time level index 
    396       REAL(wp)                        , INTENT(in   ) ::   p2dt            ! tracer time-step 
    397       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pbef            ! before field 
    398       REAL(wp), DIMENSION(A2D(nn_hls)    ,jpk), INTENT(in   ) ::   paft            ! after field 
    399       REAL(dp), DIMENSION(A2D(nn_hls)    ,jpk), INTENT(inout) ::   paa, pbb, pcc   ! monotonic fluxes in the 3 directions 
    400       ! 
    401       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    402       INTEGER  ::   ikm1         ! local integer 
    403       REAL(dp) ::   zpos, zneg, zbt, za, zb, zc, zbig, zrtrn    ! local scalars 
    404       REAL(dp) ::   zau, zbu, zcu, zav, zbv, zcv, zup, zdo            !   -      - 
    405       REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo, zbup, zbdo 
    406       !!---------------------------------------------------------------------- 
    407       ! 
    408       zbig  = 1.e+40_dp 
    409       zrtrn = 1.e-15_dp 
    410       zbetup(:,:,:) = 0._dp   ;   zbetdo(:,:,:) = 0._dp 
    411  
    412       ! Search local extrema 
    413       ! -------------------- 
    414       ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 
    415       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 
    416          zbup(ji,jj,jk) = MAX( pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) ),   & 
    417             &                  paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1._wp - tmask(ji,jj,jk) )  ) 
    418          zbdo(ji,jj,jk) = MIN( pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) ),   & 
    419             &                  paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1._wp - tmask(ji,jj,jk) )  ) 
    420       END_3D 
    421  
    422       DO jk = 1, jpkm1 
    423          ikm1 = MAX(jk-1,1) 
    424          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    425  
    426             ! search maximum in neighbourhood 
    427             zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 
    428                &        zbup(ji-1,jj  ,jk  ), zbup(ji+1,jj  ,jk  ),   & 
    429                &        zbup(ji  ,jj-1,jk  ), zbup(ji  ,jj+1,jk  ),   & 
    430                &        zbup(ji  ,jj  ,ikm1), zbup(ji  ,jj  ,jk+1)  ) 
    431  
    432             ! search minimum in neighbourhood 
    433             zdo = MIN(  zbdo(ji  ,jj  ,jk  ),   & 
    434                &        zbdo(ji-1,jj  ,jk  ), zbdo(ji+1,jj  ,jk  ),   & 
    435                &        zbdo(ji  ,jj-1,jk  ), zbdo(ji  ,jj+1,jk  ),   & 
    436                &        zbdo(ji  ,jj  ,ikm1), zbdo(ji  ,jj  ,jk+1)  ) 
    437  
    438             ! positive part of the flux 
    439             zpos = MAX( 0., paa(ji-1,jj  ,jk  ) ) - MIN( 0., paa(ji  ,jj  ,jk  ) )   & 
    440                & + MAX( 0., pbb(ji  ,jj-1,jk  ) ) - MIN( 0., pbb(ji  ,jj  ,jk  ) )   & 
    441                & + MAX( 0., pcc(ji  ,jj  ,jk+1) ) - MIN( 0., pcc(ji  ,jj  ,jk  ) ) 
    442  
    443             ! negative part of the flux 
    444             zneg = MAX( 0., paa(ji  ,jj  ,jk  ) ) - MIN( 0., paa(ji-1,jj  ,jk  ) )   & 
    445                & + MAX( 0., pbb(ji  ,jj  ,jk  ) ) - MIN( 0., pbb(ji  ,jj-1,jk  ) )   & 
    446                & + MAX( 0., pcc(ji  ,jj  ,jk  ) ) - MIN( 0., pcc(ji  ,jj  ,jk+1) ) 
    447  
    448             ! up & down beta terms 
    449             zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt 
    450             zbetup(ji,jj,jk) = ( zup            - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 
    451             zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo            ) / ( zneg + zrtrn ) * zbt 
    452          END_2D 
    453       END DO 
    454       IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp )   ! lateral boundary cond. (unchanged sign) 
    455  
    456       ! 3. monotonic flux in the i & j direction (paa & pbb) 
    457       ! ---------------------------------------- 
    458       DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
    459          zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 
    460          zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 
    461          zcu =       ( 0.5  + SIGN( 0.5_wp , CASTWP(paa(ji,jj,jk)) ) ) 
    462          paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 
    463  
    464          zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 
    465          zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 
    466          zcv =       ( 0.5  + SIGN( 0.5_wp , CASTWP(pbb(ji,jj,jk)) ) ) 
    467          pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 
    468  
    469       ! monotonic flux in the k direction, i.e. pcc 
    470       ! ------------------------------------------- 
    471          za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 
    472          zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 
    473          zc =       ( 0.5  + SIGN( 0.5_wp , CASTWP(pcc(ji,jj,jk+1)) ) ) 
    474          pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 
    475       END_3D 
    476       ! 
    477    END SUBROUTINE nonosc 
    478  
    479  
    480    SUBROUTINE interp_4th_cpt_org( pt_in, pt_out ) 
    481       !!---------------------------------------------------------------------- 
    482       !!                  ***  ROUTINE interp_4th_cpt_org  *** 
    483       !! 
    484       !! **  Purpose :   Compute the interpolation of tracer at w-point 
    485       !! 
    486       !! **  Method  :   4th order compact interpolation 
    487       !!---------------------------------------------------------------------- 
    488       REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! now tracer fields 
    489       REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   pt_out   ! now tracer field interpolated at w-pts 
    490       ! 
    491       INTEGER :: ji, jj, jk   ! dummy loop integers 
    492       REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 
    493       !!---------------------------------------------------------------------- 
    494  
    495       DO_3D( 1, 1, 1, 1, 3, jpkm1 )       !==  build the three diagonal matrix  ==! 
    496          zwd (ji,jj,jk) = 4._wp 
    497          zwi (ji,jj,jk) = 1._wp 
    498          zws (ji,jj,jk) = 1._wp 
    499          zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
    500          ! 
    501          IF( tmask(ji,jj,jk+1) == 0._wp) THEN   ! Switch to second order centered at bottom 
    502             zwd (ji,jj,jk) = 1._wp 
    503             zwi (ji,jj,jk) = 0._wp 
    504             zws (ji,jj,jk) = 0._wp 
    505             zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
    506          ENDIF 
    507       END_3D 
    508       ! 
    509       jk = 2                                    ! Switch to second order centered at top 
    510       DO_2D( 1, 1, 1, 1 ) 
    511          zwd (ji,jj,jk) = 1._wp 
    512          zwi (ji,jj,jk) = 0._wp 
    513          zws (ji,jj,jk) = 0._wp 
    514          zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 
    515       END_2D 
    516       ! 
    517       !                       !==  tridiagonal solve  ==! 
    518       DO_2D( 1, 1, 1, 1 )           ! first recurrence 
    519          zwt(ji,jj,2) = zwd(ji,jj,2) 
    520       END_2D 
    521       DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 
    522          zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
    523       END_3D 
    524       ! 
    525       DO_2D( 1, 1, 1, 1 )           ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    526          pt_out(ji,jj,2) = zwrm(ji,jj,2) 
    527       END_2D 
    528       DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 
    529          pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 
    530       END_3D 
    531  
    532       DO_2D( 1, 1, 1, 1 )           ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 
    533          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
    534       END_2D 
    535       DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 ) 
    536          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    537       END_3D 
    538       ! 
    539    END SUBROUTINE interp_4th_cpt_org 
    540  
    541  
    542    SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 
    543       !!---------------------------------------------------------------------- 
    544       !!                  ***  ROUTINE interp_4th_cpt  *** 
    545       !! 
    546       !! **  Purpose :   Compute the interpolation of tracer at w-point 
    547       !! 
    548       !! **  Method  :   4th order compact interpolation 
    549       !!---------------------------------------------------------------------- 
    550       REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pt_in    ! field at t-point 
    551       REAL(wp),DIMENSION(A2D(nn_hls)    ,jpk), INTENT(  out) ::   pt_out   ! field interpolated at w-point 
    552       ! 
    553       INTEGER ::   ji, jj, jk   ! dummy loop integers 
    554       INTEGER ::   ikt, ikb     ! local integers 
    555       REAL(wp),DIMENSION(A2D(nn_hls),jpk) :: zwd, zwi, zws, zwrm, zwt 
    556       !!---------------------------------------------------------------------- 
    557       ! 
    558       !                      !==  build the three diagonal matrix & the RHS  ==! 
    559       ! 
    560       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 )    ! interior (from jk=3 to jpk-1) 
    561          zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp                 !       diagonal 
    562          zwi (ji,jj,jk) =         wmask(ji,jj,jk)                         ! lower diagonal 
    563          zws (ji,jj,jk) =         wmask(ji,jj,jk)                         ! upper diagonal 
    564          zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk)                     &   ! RHS 
    565             &           *       ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 
    566       END_3D 
    567       ! 
    568 !!gm 
    569 !      SELECT CASE( kbc )               !* boundary condition 
    570 !      CASE( np_NH   )   ! Neumann homogeneous at top & bottom 
    571 !      CASE( np_CEN2 )   ! 2nd order centered  at top & bottom 
    572 !      END SELECT 
    573 !!gm 
    574       ! 
    575       IF ( ln_isfcav ) THEN            ! set level two values which may not be set in ISF case 
    576          zwd(:,:,2) = 1._wp  ;  zwi(:,:,2) = 0._wp  ;  zws(:,:,2) = 0._wp  ;  zwrm(:,:,2) = 0._wp 
    577       END IF 
    578       ! 
    579       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )              ! 2nd order centered at top & bottom 
    580          ikt = mikt(ji,jj) + 1            ! w-point below the 1st  wet point 
    581          ikb = MAX(mbkt(ji,jj), 2)        !     -   above the last wet point 
    582          ! 
    583          zwd (ji,jj,ikt) = 1._wp          ! top 
    584          zwi (ji,jj,ikt) = 0._wp 
    585          zws (ji,jj,ikt) = 0._wp 
    586          zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) ) 
    587          ! 
    588          zwd (ji,jj,ikb) = 1._wp          ! bottom 
    589          zwi (ji,jj,ikb) = 0._wp 
    590          zws (ji,jj,ikb) = 0._wp 
    591          zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 
    592       END_2D 
    593       ! 
    594       !                       !==  tridiagonal solver  ==! 
    595       ! 
    596       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )           !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1 
    597          zwt(ji,jj,2) = zwd(ji,jj,2) 
    598       END_2D 
    599       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 
    600          zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
    601       END_3D 
    602       ! 
    603       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    604          pt_out(ji,jj,2) = zwrm(ji,jj,2) 
    605       END_2D 
    606       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 3, jpkm1 ) 
    607          pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 
    608       END_3D 
    609  
    610       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk 
    611          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
    612       END_2D 
    613       DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 ) 
    614          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    615       END_3D 
    616       ! 
    617    END SUBROUTINE interp_4th_cpt 
    618  
    619  
    620    SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev ) 
    621       !!---------------------------------------------------------------------- 
    622       !!                  ***  ROUTINE tridia_solver  *** 
    623       !! 
    624       !! **  Purpose :   solve a symmetric 3diagonal system 
    625       !! 
    626       !! **  Method  :   solve M.t_out = RHS(t)  where M is a tri diagonal matrix ( jpk*jpk ) 
    627       !! 
    628       !!             ( D_1 U_1  0   0   0  )( t_1 )   ( RHS_1 ) 
    629       !!             ( L_2 D_2 U_2  0   0  )( t_2 )   ( RHS_2 ) 
    630       !!             (  0  L_3 D_3 U_3  0  )( t_3 ) = ( RHS_3 ) 
    631       !!             (        ...          )( ... )   ( ...  ) 
    632       !!             (  0   0   0  L_k D_k )( t_k )   ( RHS_k ) 
    633       !! 
    634       !!        M is decomposed in the product of an upper and lower triangular matrix. 
    635       !!        The tri-diagonals matrix is given as input 3D arrays:   pD, pU, pL 
    636       !!        (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal). 
    637       !!        The solution is pta. 
    638       !!        The 3d array zwt is used as a work space array. 
    639       !!---------------------------------------------------------------------- 
    640       REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pD, pU, pL    ! 3-diagonal matrix 
    641       REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   pRHS          ! Right-Hand-Side 
    642       REAL(wp),DIMENSION(A2D(nn_hls),jpk), INTENT(  out) ::   pt_out        !!gm field at level=F(klev) 
    643       INTEGER                    , INTENT(in   ) ::   klev          ! =1 pt_out at w-level 
    644       !                                                             ! =0 pt at t-level 
    645       INTEGER ::   ji, jj, jk   ! dummy loop integers 
    646       INTEGER ::   kstart       ! local indices 
    647       REAL(wp),DIMENSION(A2D(nn_hls),jpk) ::   zwt   ! 3D work array 
    648       !!---------------------------------------------------------------------- 
    649       ! 
    650       kstart =  1  + klev 
    651       ! 
    652       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                         !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1 
    653          zwt(ji,jj,kstart) = pD(ji,jj,kstart) 
    654       END_2D 
    655       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 
    656          zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 
    657       END_3D 
    658       ! 
    659       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                        !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1 
    660          pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 
    661       END_2D 
    662       DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, kstart+1, jpkm1 ) 
    663          pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 
    664       END_3D 
    665  
    666       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                       !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk 
    667          pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 
    668       END_2D 
    669       DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, kstart, -1 ) 
    670          pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 
    671       END_3D 
    672       ! 
    673    END SUBROUTINE tridia_solver 
    674  
     1008   END SUBROUTINE tra_adv_fct_lf 
     1009#endif 
    6751010   !!====================================================================== 
    6761011END MODULE traadv_fct 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traadv_mus.F90

    r14644 r14986  
    8282      LOGICAL                                  , INTENT(in   ) ::   ld_msc_ups      ! use upstream scheme within muscl 
    8383      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
    84       ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 
     84      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    8585      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume flux components 
    8686      REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    … …  
    9494      !!---------------------------------------------------------------------- 
    9595      ! 
    96       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     96      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    9797         IF( kt == kit000 )  THEN 
    9898            IF(lwp) WRITE(numout,*) 
    … …  
    140140            zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 
    141141         END_3D 
    142          ! lateral boundary conditions   (changed sign) 
    143          IF ( nn_hls.EQ.1 ) CALL lbc_lnk( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) 
    144142         !                                !-- Slopes of tracer 
    145143         zslpx(:,:,jpk) = 0._wp                 ! bottom values 
    146144         zslpy(:,:,jpk) = 0._wp 
    147          DO_3D( nn_hls-1, 1, nn_hls-1, 1, 1, jpkm1 ) 
     145         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
    148146            zslpx(ji,jj,jk) =                       ( zwx(ji,jj,jk) + zwx(ji-1,jj  ,jk) )   & 
    149147               &            * ( 0.25 + SIGN( 0.25_wp, zwx(ji,jj,jk) * zwx(ji-1,jj  ,jk) ) ) 
    … …  
    152150         END_3D 
    153151         ! 
    154          DO_3D( nn_hls-1, 1, nn_hls-1, 1, 1, jpkm1 )    !-- Slopes limitation 
     152         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )    !-- Slopes limitation 
    155153            zslpx(ji,jj,jk) = SIGN( 1.0_wp, zslpx(ji,jj,jk) ) * MIN(    ABS( zslpx(ji  ,jj,jk) ),   & 
    156154               &                                                     2.*ABS( zwx  (ji-1,jj,jk) ),   & 
    … …  
    160158               &                                                     2.*ABS( zwy  (ji,jj  ,jk) ) ) 
    161159         END_3D 
    162          ! 
    163          DO_3D( nn_hls-1, 0, nn_hls-1, 0, 1, jpkm1 )    !-- MUSCL horizontal advective fluxes 
     160         ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility 
     161         IF ( nn_hls==1 ) CALL lbc_lnk( 'traadv_mus', zslpx, 'T', -1.0_wp , zslpy, 'T', -1.0_wp )   ! lateral boundary conditions   (changed sign) 
     162         ! 
     163         DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !-- MUSCL horizontal advective fluxes 
    164164            ! MUSCL fluxes 
    165165            z0u = SIGN( 0.5_wp, pU(ji,jj,jk) ) 
    … …  
    177177            zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 
    178178         END_3D 
    179          IF ( nn_hls.EQ.1 ) CALL lbc_lnk( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )   ! lateral boundary conditions   (changed sign) 
    180179         ! 
    181180         DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !-- Tracer advective trend 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traadv_qck.F90

    r14644 r14986  
    2727   USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
    2828   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     29#if defined key_loop_fusion 
     30   USE traadv_qck_lf   ! QCK    scheme            (tra_adv_qck  routine - loop fusion version) 
     31#endif 
    2932 
    3033   IMPLICIT NONE 
    … …  
    9295      INTEGER                                  , INTENT(in   ) ::   kjpt            ! number of tracers 
    9396      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
    94       ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 
     97      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    9598      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components 
    9699      REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    97100      !!---------------------------------------------------------------------- 
    98101      ! 
    99       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     102#if defined key_loop_fusion 
     103      CALL tra_adv_qck_lf ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs ) 
     104#else 
     105      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    100106         IF( kt == kit000 )  THEN 
    101107            IF(lwp) WRITE(numout,*) 
    … …  
    118124      CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs ) 
    119125      ! 
     126#endif 
    120127   END SUBROUTINE tra_adv_qck 
    121128 
    … …  
    130137      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers 
    131138      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step 
    132       ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 
     139      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    133140      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU        ! i-velocity components 
    134141      REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    … …  
    150157            zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb)        ! Downstream in the x-direction for the tracer 
    151158         END_3D 
    152          IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions 
     159         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )   ! Lateral boundary conditions 
    153160 
    154161         ! 
    … …  
    168175         END_3D 
    169176         !--- Lateral boundary conditions 
    170          IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp ) 
     177         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp,  zwx(:,:,:), 'T', 1.0_wp ) 
    171178 
    172179         !--- QUICKEST scheme 
    … …  
    177184            zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 
    178185         END_3D 
    179          IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )      ! Lateral boundary conditions 
     186         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )      ! Lateral boundary conditions 
    180187 
    181188         ! 
    … …  
    215222      INTEGER                                  , INTENT(in   ) ::   kjpt       ! number of tracers 
    216223      REAL(wp)                                 , INTENT(in   ) ::   p2dt       ! tracer time-step 
    217       ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 
     224      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    218225      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pV        ! j-velocity components 
    219226      REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    … …  
    230237         zfd(:,:,:) = 0.0     ;   zwy(:,:,:) = 0.0 
    231238         ! 
    232          DO jk = 1, jpkm1 
    233             ! 
    234             !--- Computation of the ustream and downstream value of the tracer and the mask 
    235             DO_2D( 0, 0, nn_hls-1, nn_hls-1 ) 
    236                ! Upstream in the x-direction for the tracer 
    237                zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb) 
    238                ! Downstream in the x-direction for the tracer 
    239                zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb) 
    240             END_2D 
    241          END DO 
    242          IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp )   ! Lateral boundary conditions 
    243  
     239         !--- Computation of the ustream and downstream value of the tracer and the mask 
     240         DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
     241            ! Upstream in the x-direction for the tracer 
     242            zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb) 
     243            ! Downstream in the x-direction for the tracer 
     244            zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb) 
     245         END_3D 
     246 
     247         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )   ! Lateral boundary conditions 
     248 
     249         ! Correct zfd on northfold after lbc_lnk; see #2640 
     250         IF( nn_hls == 1 .AND. l_IdoNFold .AND. ntej == Nje0 ) THEN 
     251            DO jk = 1, jpkm1 
     252               WHERE( tmask_i(ntsi:ntei,ntej:jpj) == 0._wp ) zfd(ntsi:ntei,ntej:jpj,jk) = zfc(ntsi:ntei,ntej:jpj,jk) 
     253            END DO 
     254         ENDIF 
    244255         ! 
    245256         ! Horizontal advective fluxes 
    … …  
    260271 
    261272         !--- Lateral boundary conditions 
    262          IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp ) 
     273         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp ) 
    263274 
    264275         !--- QUICKEST scheme 
    … …  
    269280            zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 
    270281         END_3D 
    271          IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp )    !--- Lateral boundary conditions 
     282         IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. )    !--- Lateral boundary conditions 
    272283         ! 
    273284         ! Tracer flux on the x-direction 
    … …  
    307318      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    308319      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers 
    309       ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 
     320      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    310321      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pW      ! vertical velocity 
    311322      REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! active tracers and RHS of tracer equation 
    … …  
    366377      !---------------------------------------------------------------------- 
    367378      ! 
    368       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 
     379      DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
    369380         zc     = puc(ji,jj,jk)                         ! Courant number 
    370381         zcurv  = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk) 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traadv_ubs.F90

    r14644 r14986  
    2626   USE lbclnk         ! ocean lateral boundary condition (or mpp link) 
    2727   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     28#if defined key_loop_fusion 
     29   USE traadv_ubs_lf  ! UBS      scheme            (tra_adv_ubs  routine - loop fusion version) 
     30#endif 
    2831 
    2932   IMPLICIT NONE 
    … …  
    9396      INTEGER                                  , INTENT(in   ) ::   kn_ubs_v        ! number of tracers 
    9497      REAL(wp)                                 , INTENT(in   ) ::   p2dt            ! tracer time-step 
    95       ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support) 
     98      ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
    9699      REAL(wp), DIMENSION(jpi,jpj,jpk         ), INTENT(in   ) ::   pU, pV, pW      ! 3 ocean volume transport components 
    97100      REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt              ! tracers and RHS of tracer equation 
    … …  
    104107      !!---------------------------------------------------------------------- 
    105108      ! 
    106       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     109#if defined key_loop_fusion 
     110      CALL tra_adv_ubs_lf    ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs, kn_ubs_v ) 
     111#else 
     112      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    107113         IF( kt == kit000 )  THEN 
    108114            IF(lwp) WRITE(numout,*) 
    … …  
    141147            ! 
    142148         END DO 
    143          IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp )   ! Lateral boundary cond. (unchanged sgn) 
     149         IF (nn_hls==1) CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp, ld4only= .TRUE. )   ! Lateral boundary cond. (unchanged sgn) 
    144150         ! 
    145151         DO_3D( 1, 0, 1, 0, 1, jpkm1 )   !==  Horizontal advective fluxes  ==!     (UBS) 
    … …  
    156162         END_3D 
    157163         ! 
    158          DO_3D( 1, 1, 1, 1, 1, jpk ) 
     164         DO_3D( 0, 0, 0, 0, 1, jpk ) 
    159165            zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs)      ! store the initial trends before its update 
    160166         END_3D 
    … …  
    170176         END DO 
    171177         ! 
    172          DO_3D( 1, 1, 1, 1, 1, jpk ) 
     178         DO_3D( 0, 0, 0, 0, 1, jpk ) 
    173179            zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltu(ji,jj,jk)  ! Horizontal advective trend used in vertical 2nd order FCT case 
    174180         END_3D                                                     ! and/or in trend diagnostic (l_trd=T) 
    … …  
    198204            ! 
    199205            !                               !*  upstream advection with initial mass fluxes & intermediate update  ==! 
    200             DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 
     206            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    201207               zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 
    202208               zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 
    … …  
    205211            IF( ln_linssh ) THEN                ! top ocean value (only in linear free surface as ztw has been w-masked) 
    206212               IF( ln_isfcav ) THEN                   ! top of the ice-shelf cavities and at the ocean surface 
    207                   DO_2D( 1, 1, 1, 1 ) 
     213                  DO_2D( 0, 0, 0, 0 ) 
    208214                     ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb)   ! linear free surface 
    209215                  END_2D 
    210216               ELSE                                   ! no cavities: only at the ocean surface 
    211                   DO_2D( 1, 1, 1, 1 ) 
     217                  DO_2D( 0, 0, 0, 0 ) 
    212218                     ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) 
    213219                  END_2D 
    … …  
    223229            ! 
    224230            !                          !*  anti-diffusive flux : high order minus low order 
    225             DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 
     231            DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    226232               ztw(ji,jj,jk) = (   0.5_wp * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) )   & 
    227233                  &              - ztw(ji,jj,jk)   ) * wmask(ji,jj,jk) 
    … …  
    238244            END_3D 
    239245            IF( ln_linssh ) THEN 
    240                DO_2D( 1, 1, 1, 1 ) 
     246               DO_2D( 0, 0, 0, 0 ) 
    241247                  ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)     !!gm ISF & 4th COMPACT doesn't work 
    242248               END_2D 
    … …  
    261267      END DO 
    262268      ! 
     269#endif 
    263270   END SUBROUTINE tra_adv_ubs 
    264271 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/trabbc.F90

    r14219 r14986  
    103103      ENDIF 
    104104      ! 
    105       IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain 
    106          CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) ) 
    107       ENDIF 
     105      CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) ) 
     106 
    108107      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=CASTWP(pts(:,:,:,jp_tem,Krhs)), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 
    109108      ! 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/trabbl.F90

    r14644 r14986  
    127127         CALL prt_ctl( tab3d_1=CASTWP(pts(:,:,:,jp_tem,Krhs)), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, & 
    128128            &          tab3d_2=CASTWP(pts(:,:,:,jp_sal,Krhs)), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    129          IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile 
    130             CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef 
    131             CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef 
    132          ENDIF 
     129         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef 
     130         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef 
    133131         ! 
    134132      ENDIF 
    … …  
    140138         CALL prt_ctl( tab3d_1=CASTWP(pts(:,:,:,jp_tem,Krhs)), clinfo1=' bbl_adv  - Ta: ', mask1=tmask, & 
    141139            &          tab3d_2=CASTWP(pts(:,:,:,jp_sal,Krhs)), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    142          IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile 
    143             ! lateral boundary conditions ; just need for outputs 
    144             CALL lbc_lnk( 'trabbl', utr_bbl, 'U', 1.0_wp , vtr_bbl, 'V', 1.0_wp ) 
    145             CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport 
    146             CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport 
    147          ENDIF 
     140         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport 
     141         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport 
    148142         ! 
    149143      ENDIF 
    … …  
    216210 
    217211 
     212   ! NOTE: [tiling] tiling changes the results, but only the order of floating point operations is different 
    218213   SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm ) 
    219214      !!---------------------------------------------------------------------- 
    … …  
    239234      INTEGER  ::   iis , iid , ijs , ijd    ! local integers 
    240235      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       - 
    241       INTEGER  ::   isi, isj                 !   -       - 
    242236      REAL(wp) ::   zbtr, ztra               ! local scalars 
    243237      REAL(wp) ::   zu_bbl, zv_bbl           !   -      - 
    244238      !!---------------------------------------------------------------------- 
    245       ! 
    246       IF( ntsi == Nis0 ) THEN ; isi = 1 ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling 
    247       IF( ntsj == Njs0 ) THEN ; isj = 1 ; ELSE ; isj = 0 ; ENDIF 
    248239      !                                                          ! =========== 
    249240      DO jn = 1, kjpt                                            ! tracer loop 
    250241         !                                                       ! =========== 
    251          DO_2D( isi, 0, isj, 0 )            ! CAUTION start from i=1 to update i=2 when cyclic east-west 
     242         DO_2D_OVR( 1, 0, 1, 0 )            ! CAUTION start from i=1 to update i=2 when cyclic east-west 
    252243            IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection 
    253244               ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf) 
    … …  
    341332      !!---------------------------------------------------------------------- 
    342333      ! 
    343       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     334      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    344335         IF( kt == kit000 )  THEN 
    345336            IF(lwp)  WRITE(numout,*) 
    … …  
    364355      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   ! 
    365356         !                                !-------------------! 
    366          DO_2D( 1, 0, 1, 0 )                   ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) 
     357         DO_2D_OVR( 1, 0, 1, 0 )                   ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) 
    367358            !                                                   ! i-direction 
    368359            za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point 
    … …  
    394385         ! 
    395386         CASE( 1 )                                   != use of upper velocity 
    396             DO_2D( 1, 0, 1, 0 )                              ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0 
     387            DO_2D_OVR( 1, 0, 1, 0 )                              ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0 
    397388               !                                                  ! i-direction 
    398389               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
    … …  
    423414         CASE( 2 )                                 != bbl velocity = F( delta rho ) 
    424415            zgbbl = grav * rn_gambbl 
    425             DO_2D( 1, 0, 1, 0 )                         ! criteria: rho_up > rho_down 
     416            DO_2D_OVR( 1, 0, 1, 0 )                         ! criteria: rho_up > rho_down 
    426417               !                                                  ! i-direction 
    427418               ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf) 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/tradmp.F90

    r14286 r14986  
    5454#  include "do_loop_substitute.h90" 
    5555#  include "single_precision_substitute.h90" 
     56#  include "domzgr_substitute.h90" 
    5657   !!---------------------------------------------------------------------- 
    5758   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    9798      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices 
    9899      REAL(dp), DIMENSION(A2D(nn_hls),jpk,jpts)     ::  zts_dta 
     100      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE ::  zwrk 
    99101      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::  ztrdts 
    100102      !!---------------------------------------------------------------------- 
    … …  
    102104      IF( ln_timing )   CALL timing_start('tra_dmp') 
    103105      ! 
    104       IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    105          ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 
    106          ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) 
     106      IF( l_trdtra .OR. iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN   !* Save ta and sa trends 
     107         ALLOCATE( ztrdts(A2D(nn_hls),jpk,jpts) ) 
     108         DO jn = 1, jpts 
     109            DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 
     110               ztrdts(ji,jj,jk,jn) = pts(ji,jj,jk,jn,Krhs) 
     111            END_3D 
     112         END DO 
    107113      ENDIF 
    108114      !                           !==  input T-S data at kt  ==! 
    … …  
    140146         ! 
    141147      END SELECT 
     148      ! 
     149      ! outputs (clem trunk) 
     150      IF( iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN 
     151         ALLOCATE( zwrk(A2D(nn_hls),jpk) )          ! Needed to handle expressions containing e3t when using key_qco or key_linssh 
     152         zwrk(:,:,:) = 0._wp 
     153 
     154         IF( iom_use('hflx_dmp_cea') ) THEN 
     155            DO_3D( 0, 0, 0, 0, 1, jpk ) 
     156               zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_tem,Krhs) - ztrdts(ji,jj,jk,jp_tem) ) * e3t(ji,jj,jk,Kmm) 
     157            END_3D 
     158            CALL iom_put('hflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rcp * rho0 ) ! W/m2 
     159         ENDIF 
     160         IF( iom_use('sflx_dmp_cea') ) THEN 
     161            DO_3D( 0, 0, 0, 0, 1, jpk ) 
     162               zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_sal,Krhs) - ztrdts(ji,jj,jk,jp_sal) ) * e3t(ji,jj,jk,Kmm) 
     163            END_3D 
     164            CALL iom_put('sflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rho0 )       ! g/m2/s 
     165         ENDIF 
     166 
     167         DEALLOCATE( zwrk ) 
     168      ENDIF 
    142169      ! 
    143170      IF( l_trdtra )   THEN       ! trend diagnostic 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traisf.F90

    r14219 r14986  
    4848      IF( ln_timing )   CALL timing_start('tra_isf') 
    4949      ! 
    50       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     50      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    5151         IF( kt == nit000 ) THEN 
    5252            IF(lwp) WRITE(numout,*) 
    … …  
    8080      ! 
    8181      IF ( ln_isfdebug ) THEN 
    82          IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only for the full domain 
     82         IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only for the full domain 
    8383            CALL debug('tra_isf: pts(:,:,:,:,Krhs) T', CASTWP(pts(:,:,:,1,Krhs))) 
    8484            CALL debug('tra_isf: pts(:,:,:,:,Krhs) S', CASTWP(pts(:,:,:,2,Krhs))) 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traldf.F90

    r14219 r14986  
    1717   USE oce            ! ocean dynamics and tracers 
    1818   USE dom_oce        ! ocean space and time domain 
    19    ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 
    20    USE domtile 
    2119   USE phycst         ! physical constants 
    2220   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff. 
    … …  
    5957      !! 
    6058      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds 
    61       ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 
    62       LOGICAL :: lskip 
    6359      !!---------------------------------------------------------------------- 
    6460      ! 
    6561      IF( ln_timing )   CALL timing_start('tra_ldf') 
    6662      ! 
    67       lskip = .FALSE. 
    68  
    6963      IF( l_trdtra )   THEN                    !* Save ta and sa trends 
    7064         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
    … …  
    7266         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 
    7367      ENDIF 
    74  
    75       ! TEMP: [tiling] These changes not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 
    76       IF( nldf_tra == np_blp .OR. nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it )  THEN 
    77          IF( ln_tile ) THEN 
    78             IF( ntile == 1 ) THEN 
    79                CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 
    80             ELSE 
    81                lskip = .TRUE. 
    82             ENDIF 
    83          ENDIF 
    84       ENDIF 
    85       IF( .NOT. lskip ) THEN 
    86          ! 
    87          SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend 
    88          CASE ( np_lap   )                                  ! laplacian: iso-level operator 
    89             CALL tra_ldf_lap  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, CASTWP(pts(:,:,:,:,Kbb)), pts(:,:,:,:,Krhs),                   jpts,  1 ) 
    90          CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec) 
    91             CALL tra_ldf_iso  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, CASTWP(pts(:,:,:,:,Kbb)), CASTWP(pts(:,:,:,:,Kbb)), pts(:,:,:,:,Krhs), jpts,  1 ) 
    92          CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies) 
    93             CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, CASTWP(pts(:,:,:,:,Kbb)), CASTWP(pts(:,:,:,:,Kbb)), pts(:,:,:,:,Krhs), jpts,  1 ) 
    94          CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators 
    95             IF(nn_hls.EQ.2) CALL lbc_lnk( 'tra_ldf', pts(:,:,:,:,Kbb), 'T',1._wp) 
    96             CALL tra_ldf_blp  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, CASTWP(pts(:,:,:,:,Kbb)), pts(:,:,:,:,Krhs),             jpts, nldf_tra ) 
    97          END SELECT 
    98          ! 
    99          IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics 
    100             ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
    101             ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 
    102             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_ldf, ztrdt ) 
    103             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds ) 
    104             DEALLOCATE( ztrdt, ztrds ) 
    105          ENDIF 
    106  
    107          ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 
    108          IF( ln_tile .AND. ntile == 0 ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) 
     68      ! 
     69      SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend 
     70      CASE ( np_lap   )                                  ! laplacian: iso-level operator 
     71         CALL tra_ldf_lap  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, CASTWP(pts(:,:,:,:,Kbb)), pts(:,:,:,:,Krhs),                   jpts,  1 ) 
     72      CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec) 
     73         CALL tra_ldf_iso  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, CASTWP(pts(:,:,:,:,Kbb)), CASTWP(pts(:,:,:,:,Kbb)), pts(:,:,:,:,Krhs), jpts,  1 ) 
     74      CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies) 
     75         CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, CASTWP(pts(:,:,:,:,Kbb)), CASTWP(pts(:,:,:,:,Kbb)), pts(:,:,:,:,Krhs), jpts,  1 ) 
     76      CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators 
     77         CALL tra_ldf_blp  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, CASTWP(pts(:,:,:,:,Kbb)), pts(:,:,:,:,Krhs),             jpts, nldf_tra ) 
     78      END SELECT 
     79      ! 
     80      IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics 
     81         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 
     82         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 
     83         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_ldf, ztrdt ) 
     84         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds ) 
     85         DEALLOCATE( ztrdt, ztrds ) 
    10986      ENDIF 
    11087      !                                        !* print mean trends (used for debugging) 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traldf_iso.F90

    r14221 r14986  
    132132      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    133133      INTEGER  ::  ikt 
    134       INTEGER  ::  ierr             ! local integer 
     134      INTEGER  ::  ierr, iij        ! local integer 
    135135      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars 
    136136      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      - 
    … …  
    141141      ! 
    142142      IF( kpass == 1 .AND. kt == kit000 )  THEN 
    143          IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     143         IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    144144            IF(lwp) WRITE(numout,*) 
    145145            IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype 
    … …  
    147147         ENDIF 
    148148         ! 
    149          DO_3D( 0, 0, 0, 0, 1, jpk ) 
     149         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
    150150            akz     (ji,jj,jk) = 0._wp 
    151151            ah_wslp2(ji,jj,jk) = 0._wp 
    … …  
    153153      ENDIF 
    154154      ! 
    155       IF( ntile == 0 .OR. ntile == 1 )  THEN                           ! Do only on the first tile 
     155      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                           ! Do only on the first tile 
    156156         l_hst = .FALSE. 
    157157         l_ptr = .FALSE. 
    … …  
    161161      ENDIF 
    162162      ! 
     163      ! Define pt_rhs halo points for multi-point haloes in bilaplacian case 
     164      IF( nldf_tra == np_blp_i .AND. kpass == 1 ) THEN ; iij = nn_hls 
     165      ELSE                                             ; iij = 1 
     166      ENDIF 
     167 
    163168      ! 
    164169      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
    … …  
    172177      IF( kpass == 1 ) THEN                  !==  first pass only  ==! 
    173178         ! 
    174          DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     179         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    175180            ! 
    176181            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
    … …  
    179184               &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
    180185               ! 
    181             zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    & 
    182                &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku 
    183             zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    & 
    184                &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv 
     186            ! round brackets added to fix the order of floating point operations 
     187            ! needed to ensure halo 1 - halo 2 compatibility 
     188            zahu_w = ( (  pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)                    & 
     189               &       )                                                           & ! bracket for halo 1 - halo 2 compatibility 
     190               &       + ( pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)                   & 
     191               &         )                                                         & ! bracket for halo 1 - halo 2 compatibility 
     192               &     ) * zmsku 
     193            zahv_w = ( (  pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)                    & 
     194               &       )                                                           & ! bracket for halo 1 - halo 2 compatibility 
     195               &       + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)                   & 
     196               &         )                                                         & ! bracket for halo 1 - halo 2 compatibility 
     197               &     ) * zmskv 
    185198               ! 
    186199            ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    … …  
    189202         ! 
    190203         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient 
    191             DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     204            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
     205               ! round brackets added to fix the order of floating point operations 
     206               ! needed to ensure halo 1 - halo 2 compatibility 
    192207               akz(ji,jj,jk) = 0.25_wp * (                                                                     & 
    193                   &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
     208                  &            ( ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   & 
    194209                  &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   & 
    195                   &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   & 
    196                   &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   ) 
     210                  &            )                                                                               & ! bracket for halo 1 - halo 2 compatibility 
     211                  &            + ( ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) ) & 
     212                  &              + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) & 
     213                  &              )                                                                             & ! bracket for halo 1 - halo 2 compatibility 
     214                  &                      ) 
    197215            END_3D 
    198216            ! 
    199217            IF( ln_traldf_blp ) THEN                ! bilaplacian operator 
    200                DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     218               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    201219                  akz(ji,jj,jk) = 16._wp   & 
    202220                     &   * ah_wslp2   (ji,jj,jk)   & 
    … …  
    206224               END_3D 
    207225            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
    208                DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     226               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    209227                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
    210228                  zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    … …  
    214232           ! 
    215233         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
    216             DO_3D( 0, 0, 0, 0, 1, jpk ) 
     234            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
    217235               akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 
    218236            END_3D 
    … …  
    227245         !!   I - masked horizontal derivative 
    228246         !!---------------------------------------------------------------------- 
    229 !!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient.... 
    230          zdit (ntsi-nn_hls,:,:) = 0._wp     ;     zdit (ntei+nn_hls,:,:) = 0._wp 
    231          zdjt (ntsi-nn_hls,:,:) = 0._wp     ;     zdjt (ntei+nn_hls,:,:) = 0._wp 
    232          !!end 
     247         zdit(:,:,:) = 0._wp 
     248         zdjt(:,:,:) = 0._wp 
    233249 
    234250         ! Horizontal tracer gradient 
    235          DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
     251         DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) 
    236252            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    237253            zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    238254         END_3D 
    239255         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient 
    240             DO_2D( 1, 0, 1, 0 )           ! bottom correction (partial bottom cell) 
     256            DO_2D( iij, iij-1, iij, iij-1 )            ! bottom correction (partial bottom cell) 
    241257               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
    242258               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    243259            END_2D 
    244260            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity 
    245                DO_2D( 1, 0, 1, 0 ) 
     261               DO_2D( iij, iij-1, iij, iij-1 ) 
    246262                  IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 
    247263                  IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 
    … …  
    256272         DO jk = 1, jpkm1                                 ! Horizontal slab 
    257273            ! 
    258             DO_2D( 1, 1, 1, 1 ) 
     274            DO_2D( iij, iij, iij, iij ) 
    259275               !                             !== Vertical tracer gradient 
    260276               zdk1t(ji,jj) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)     ! level jk+1 
    … …  
    265281            END_2D 
    266282            ! 
    267             DO_2D( 1, 0, 1, 0 )           !==  Horizontal fluxes 
     283            DO_2D( iij, iij-1, iij, iij-1 )           !==  Horizontal fluxes 
    268284               zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) 
    269285               zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 
    … …  
    278294               zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 
    279295               ! 
    280                zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
    281                   &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
    282                   &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
    283                zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
    284                   &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
    285                   &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk) 
     296               ! round brackets added to fix the order of floating point operations 
     297               ! needed to ensure halo 1 - halo 2 compatibility 
     298               zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)                       & 
     299                  &               + zcof1 * ( ( zdkt (ji+1,jj) + zdk1t(ji,jj)    & 
     300                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility 
     301                  &                         + ( zdk1t(ji+1,jj) + zdkt (ji,jj)    & 
     302                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility 
     303                  &                         ) ) * umask(ji,jj,jk) 
     304               zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)                        & 
     305                  &              + zcof2 * ( ( zdkt (ji,jj+1) + zdk1t(ji,jj)     & 
     306                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility 
     307                  &                         + ( zdk1t(ji,jj+1) + zdkt (ji,jj)    & 
     308                  &                           )                                  & ! bracket for halo 1 - halo 2 compatibility 
     309                  &                         ) ) * vmask(ji,jj,jk) 
    286310            END_2D 
    287311            ! 
    288             DO_2D( 0, 0, 0, 0 )           !== horizontal divergence and add to pta 
    289                pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    & 
    290                   &       + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   & 
    291                   &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
     312            DO_2D( iij-1, iij-1, iij-1, iij-1 )           !== horizontal divergence and add to pta 
     313               ! round brackets added to fix the order of floating point operations 
     314               ! needed to ensure halo 1 - halo 2 compatibility 
     315               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)                         & 
     316                  &       + zsign * ( ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk)        & 
     317                  &                   )                                          & ! bracket for halo 1 - halo 2 compatibility 
     318                  &                 + ( zftv(ji,jj,jk) - zftv(ji,jj-1,jk)        & 
     319                  &                   )                                          & ! bracket for halo 1 - halo 2 compatibility 
     320                  &                 ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 
    292321            END_2D 
    293322         END DO                                        !   End of slab 
    … …  
    302331         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp 
    303332 
    304          DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior (2=<jk=<jpk-1) 
     333         DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 )    ! interior (2=<jk=<jpk-1) 
    305334            ! 
    306335            zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
    … …  
    317346            zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 
    318347            ! 
    319             ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      & 
    320                &                        + zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)  )   & 
    321                &           + zcoef4 * (   zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)      & 
    322                &                        + zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)  ) 
     348            ! round brackets added to fix the order of floating point operations 
     349            ! needed to ensure halo 1 - halo 2 compatibility 
     350            ztfw(ji,jj,jk) = zcoef3 * ( ( zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)    & 
     351                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility 
     352                  &                   + ( zdit(ji-1,jj  ,jk-1) + zdit(ji  ,jj  ,jk)    & 
     353                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility 
     354                  &                   )                                                & 
     355                  &        + zcoef4 * ( ( zdjt(ji  ,jj  ,jk-1) + zdjt(ji  ,jj-1,jk)    & 
     356                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility 
     357                  &                   + ( zdjt(ji  ,jj-1,jk-1) + zdjt(ji  ,jj  ,jk)    & 
     358                  &                     )                                              & ! bracket for halo 1 - halo 2 compatibility 
     359                  &                   ) 
    323360         END_3D 
    324361         !                                !==  add the vertical 33 flux  ==! 
    325362         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
    326             DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     363            DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 
    327364               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   & 
    328365                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )               & 
    … …  
    333370            SELECT CASE( kpass ) 
    334371            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    335                DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     372               DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 
    336373                  ztfw(ji,jj,jk) =   & 
    337374                     &  ztfw(ji,jj,jk) + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
    … …  
    347384         ENDIF 
    348385         ! 
    349          DO_3D( 0, 0, 0, 0, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==! 
     386         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )    !==  Divergence of vertical fluxes added to pta  ==! 
    350387            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * r1_e1e2t(ji,jj)   & 
    351388               &                                             / e3t(ji,jj,jk,Kmm) 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traldf_lap_blp.F90

    r14644 r14986  
    104104      ! 
    105105      INTEGER  ::   ji, jj, jk, jn      ! dummy loop indices 
    106       INTEGER  ::   isi, iei, isj, iej  ! local integers 
     106      INTEGER  ::   iij 
    107107      REAL(wp) ::   zsign               ! local scalars 
    108108      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   ztu, ztv, zaheeu, zaheev 
    109109      !!---------------------------------------------------------------------- 
    110110      ! 
    111       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     111      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    112112         IF( kt == nit000 .AND. lwp )  THEN 
    113113            WRITE(numout,*) 
    … …  
    123123      ENDIF 
    124124      ! 
    125       l_hst = .FALSE. 
    126       l_ptr = .FALSE. 
    127       IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtldf' ) .OR. iom_use( 'sopstldf' ) ) )     l_ptr = .TRUE. 
    128       IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 
    129          &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )  l_hst = .TRUE. 
    130       ! 
     125      ! Define pt_rhs halo points for multi-point haloes in bilaplacian case 
     126      IF( nldf_tra == np_blp .AND. kpass == 1 ) THEN ; iij = nn_hls 
     127      ELSE                                           ; iij = 1 
     128      ENDIF 
     129 
    131130      !                                !==  Initialization of metric arrays used for all tracers  ==! 
    132131      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
    … …  
    134133      ENDIF 
    135134 
    136       IF( ntsi == Nis0 ) THEN ; isi = nn_hls - 1 ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling 
    137       IF( ntsj == Njs0 ) THEN ; isj = nn_hls - 1 ; ELSE ; isj = 0 ; ENDIF 
    138       IF( ntei == Nie0 ) THEN ; iei = nn_hls - 1 ; ELSE ; iei = 0 ; ENDIF 
    139       IF( ntej == Nje0 ) THEN ; iej = nn_hls - 1 ; ELSE ; iej = 0 ; ENDIF 
    140  
    141       DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )            !== First derivative (gradient)  ==! 
     135      DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 )            !== First derivative (gradient)  ==! 
    142136         zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)   !!gm   * umask(ji,jj,jk) pah masked! 
    143137         zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)   !!gm   * vmask(ji,jj,jk) 
    … …  
    148142         !                          ! =========== ! 
    149143         ! 
    150          DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )            !== First derivative (gradient)  ==! 
     144         DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 )            !== First derivative (gradient)  ==! 
    151145            ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) 
    152146            ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) 
    153147         END_3D 
    154148         IF( ln_zps ) THEN                             ! set gradient at bottom/top ocean level 
    155             DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )                              ! bottom 
     149            DO_2D( iij, iij-1, iij, iij-1 )                              ! bottom 
    156150               ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn) 
    157151               ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn) 
    158152            END_2D 
    159153            IF( ln_isfcav ) THEN                             ! top in ocean cavities only 
    160                DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
     154               DO_2D( iij, iij-1, iij, iij-1 ) 
    161155                  IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 
    162156                  IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 
    … …  
    165159         ENDIF 
    166160         ! 
    167          DO_3D( isi, iei, isj, iej, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==! 
    168             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     & 
    169                &                                      +    ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   & 
    170                &                                      / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 
     161         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )            !== Second derivative (divergence) added to the general tracer trends  ==! 
     162            ! round brackets added to fix the order of floating point operations 
     163            ! needed to ensure halo 1 - halo 2 compatibility 
     164            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ( ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)    & 
     165               &                                          )                                    & ! bracket for halo 1 - halo 2 compatibility 
     166               &                                      +   ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk)    & 
     167               &                                          )                                    & ! bracket for halo 1 - halo 2 compatibility 
     168               &                                        ) / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 
    171169         END_3D 
    172170         ! 
    … …  
    218216      !!--------------------------------------------------------------------- 
    219217      ! 
    220       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     218      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    221219         IF( kt == kit000 .AND. lwp )  THEN 
    222220            WRITE(numout,*) 
    … …  
    242240      END SELECT 
    243241      ! 
    244       CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp )     ! Lateral boundary conditions (unchanged sign) 
     242      IF (nn_hls==1) CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp )     ! Lateral boundary conditions (unchanged sign) 
    245243      !                                               ! Partial top/bottom cell: GRADh( zlap ) 
    246244      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traldf_triad.F90

    r14644 r14986  
    1313   USE oce            ! ocean dynamics and active tracers 
    1414   USE dom_oce        ! ocean space and time domain 
    15    ! TEMP: [tiling] This change not necessary if XIOS has subdomain support 
    16    USE domtile 
    1715   USE domutl, ONLY : is_tile 
    1816   USE phycst         ! physical constants 
    … …  
    109107      REAL(dp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) ::   pt_rhs     ! tracer trend 
    110108      ! 
    111       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    112       INTEGER  ::  ip,jp,kp         ! dummy loop indices 
    113       INTEGER  ::  ierr            ! local integer 
    114       REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3    ! local scalars 
    115       REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4    !   -      - 
     109      INTEGER  ::  ji, jj, jk, jn, kp, iij   ! dummy loop indices 
    116110      REAL(wp) ::  zcoef0, ze3w_2, zsign          !   -      - 
    117111      ! 
    118       REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv 
    119       REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt 
    120       REAL(wp) ::   zah, zah_slp, zaei_slp 
    121       REAL(wp), DIMENSION(A2D(nn_hls),0:1)     ::   zdkt3d                         ! vertical tracer gradient at 2 levels 
    122       REAL(wp), DIMENSION(A2D(nn_hls)        ) ::   z2d                            ! 2D workspace 
    123       REAL(wp), DIMENSION(A2D(nn_hls)    ,jpk) ::   zdit, zdjt, zftu, zftv, ztfw   ! 3D     - 
    124       ! TEMP: [tiling] This can be A2D(nn_hls) if XIOS has subdomain support 
    125       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpsi_uw, zpsi_vw 
     112      REAL(wp) ::   zslope2, zbu, zbv, zbu1, zbv1, zslope21, zah, zah1, zah_ip1, zah_jp1, zbu_ip1, zbv_jp1 
     113      REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt, zdyt_jp1, ze3wr_jp1, zdzt_jp1, zah_slp1, zah_slp_jp1, zaei_slp_jp1 
     114      REAL(wp) ::   zah_slp, zaei_slp, zdxt_ip1, ze3wr_ip1, zdzt_ip1, zah_slp_ip1, zaei_slp_ip1, zaei_slp1 
     115      REAL(wp), DIMENSION(A2D(nn_hls),0:1) ::   zdkt3d                                           ! vertical tracer gradient at 2 levels 
     116      REAL(wp), DIMENSION(A2D(nn_hls)    ) ::   z2d                                              ! 2D workspace 
     117      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     - 
    126118      !!---------------------------------------------------------------------- 
    127119      ! 
    128       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     120      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    129121         IF( kpass == 1 .AND. kt == kit000 )  THEN 
    130122            IF(lwp) WRITE(numout,*) 
    … …  
    142134      ENDIF 
    143135      ! 
     136      ! Define pt_rhs halo points for multi-point haloes in bilaplacian case 
     137      IF( nldf_tra == np_blp_it .AND. kpass == 1 ) THEN ; iij = nn_hls 
     138      ELSE                                              ; iij = 1 
     139      ENDIF 
     140 
     141      ! 
    144142      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
    145143      ELSE                    ;   zsign = -1._wp 
    … …  
    152150      IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==! 
    153151         ! 
    154          DO_3D( 0, 0, 0, 0, 1, jpk ) 
     152         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
    155153            akz     (ji,jj,jk) = 0._wp 
    156154            ah_wslp2(ji,jj,jk) = 0._wp 
    157155         END_3D 
    158156         ! 
    159          DO ip = 0, 1                            ! i-k triads 
    160             DO kp = 0, 1 
    161                DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    162                   ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 
    163                   zbu   = e1e2u(ji-ip,jj) * e3u(ji-ip,jj,jk,Kmm) 
    164                   zah   = 0.25_wp * pahu(ji-ip,jj,jk) 
    165                   zslope_skew = triadi_g(ji,jj,jk,1-ip,kp) 
    166                   ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 
    167                   zslope2 = zslope_skew + ( gdept(ji-ip+1,jj,jk,Kmm) - gdept(ji-ip,jj,jk,Kmm) ) * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 
    168                   zslope2 = zslope2 *zslope2 
    169                   ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 
    170                   akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + zah * r1_e1u(ji-ip,jj)       & 
    171                      &                                                      * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 
    172                      ! 
    173                END_3D 
    174             END DO 
     157         DO kp = 0, 1                            ! i-k triads 
     158            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
     159               ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 
     160               zbu   = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 
     161               zbu1  = e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) 
     162               zah   = 0.25_wp * pahu(ji,jj,jk) 
     163               zah1  = 0.25_wp * pahu(ji-1,jj,jk) 
     164               ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 
     165               zslope2 = triadi_g(ji,jj,jk,1,kp) + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 
     166               zslope2 = zslope2 *zslope2 
     167               zslope21 = triadi_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji-1,jj,jk,Kmm) ) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp) 
     168               zslope21 = zslope21 *zslope21 
     169               ! round brackets added to fix the order of floating point operations 
     170               ! needed to ensure halo 1 - halo 2 compatibility 
     171               ah_wslp2(ji,jj,jk+kp) =  ah_wslp2(ji,jj,jk+kp) + ( zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2                    & 
     172                        &                                       + zah1 * zbu1 * ze3wr * r1_e1e2t(ji,jj) * zslope21                 & 
     173                        &                                       )                                                                  ! bracket for halo 1 - halo 2 compatibility 
     174               akz     (ji,jj,jk+kp) =  akz     (ji,jj,jk+kp) + ( zah * r1_e1u(ji,jj) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)         & 
     175                                                                + zah1 * r1_e1u(ji-1,jj) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp)  & 
     176                        &                                       )                                                                  ! bracket for halo 1 - halo 2 compatibility 
     177            END_3D 
    175178         END DO 
    176179         ! 
    177          DO jp = 0, 1                            ! j-k triads 
    178             DO kp = 0, 1 
    179                DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    180                   ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm) 
    181                   zbv   = e1e2v(ji,jj-jp) * e3v(ji,jj-jp,jk,Kmm) 
    182                   zah   = 0.25_wp * pahv(ji,jj-jp,jk) 
    183                   zslope_skew = triadj_g(ji,jj,jk,1-jp,kp) 
    184                   ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 
    185                   !    (do this by *adding* gradient of depth) 
    186                   zslope2 = zslope_skew + ( gdept(ji,jj-jp+1,jk,Kmm) - gdept(ji,jj-jp,jk,Kmm) ) * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 
    187                   zslope2 = zslope2 * zslope2 
    188                   ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 
    189                   akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + zah * r1_e2v(ji,jj-jp)     & 
    190                      &                                                      * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 
    191                   ! 
    192                END_3D 
    193             END DO 
     180         DO kp = 0, 1                            ! j-k triads 
     181            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 
     182               ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm) 
     183               zbv   = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     184               zbv1   = e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) 
     185               zah   = 0.25_wp * pahv(ji,jj,jk) 
     186               zah1   = 0.25_wp * pahv(ji,jj-1,jk) 
     187               ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 
     188               !    (do this by *adding* gradient of depth) 
     189               zslope2 = triadj_g(ji,jj,jk,1,kp) + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 
     190               zslope2 = zslope2 * zslope2 
     191               zslope21 = triadj_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji,jj-1,jk,Kmm) ) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp) 
     192               zslope21 = zslope21 * zslope21 
     193               ! round brackets added to fix the order of floating point operations 
     194               ! needed to ensure halo 1 - halo 2 compatibility 
     195               ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2                     & 
     196                        &                                      + zah1 * zbv1 * ze3wr * r1_e1e2t(ji,jj) * zslope21                  & 
     197                        &                                      )                                                                   ! bracket for halo 1 - halo 2 compatibility 
     198               akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + ( zah * r1_e2v(ji,jj) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)          & 
     199                        &                                      + zah1 * r1_e2v(ji,jj-1) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp)   & 
     200                        &                                      )                                                                   ! bracket for halo 1 - halo 2 compatibility 
     201            END_3D 
    194202         END DO 
    195203         ! 
    … …  
    197205            ! 
    198206            IF( ln_traldf_blp ) THEN                ! bilaplacian operator 
    199                DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     207               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    200208                  akz(ji,jj,jk) = 16._wp           & 
    201209                     &   * ah_wslp2   (ji,jj,jk)   & 
    … …  
    205213               END_3D 
    206214            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
    207                DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     215               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    208216                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
    209217                  zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
    … …  
    213221           ! 
    214222         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
    215             DO_3D( 0, 0, 0, 0, 1, jpk ) 
     223            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
    216224               akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 
    217225            END_3D 
    218226         ENDIF 
    219227         ! 
    220          ! TEMP: [tiling] These changes not necessary if XIOS has subdomain support 
    221          IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain 
    222             IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 
    223                IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 
    224  
    225                zpsi_uw(:,:,:) = 0._wp 
    226                zpsi_vw(:,:,:) = 0._wp 
    227  
    228                DO jp = 0, 1 
    229                   DO kp = 0, 1 
    230                      DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
    231                         zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 
    232                            & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp) 
    233                         zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 
    234                            & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+jp,jk,1-jp,kp) 
    235                      END_3D 
    236                   END DO 
    237                END DO 
    238                CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 
    239  
    240                IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile ) 
    241             ENDIF 
     228         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 
     229            zpsi_uw(:,:,:) = 0._wp 
     230            zpsi_vw(:,:,:) = 0._wp 
     231 
     232            DO kp = 0, 1 
     233               DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
     234                  ! round brackets added to fix the order of floating point operations 
     235                  ! needed to ensure halo 1 - halo 2 compatibility 
     236                  zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)                                     & 
     237                     & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji,jj,jk,1,kp)        & 
     238                     &   + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp)      & 
     239                     &   )                                                                        ! bracket for halo 1 - halo 2 compatibility 
     240                  zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp)                                     & 
     241                     & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj,jk,1,kp)        & 
     242                     &   + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp)      & 
     243                     &   )                                                                        ! bracket for halo 1 - halo 2 compatibility 
     244               END_3D 
     245            END DO 
     246            CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 
    242247         ENDIF 
    243248         ! 
    … …  
    252257         zftu(:,:,:) = 0._wp 
    253258         zftv(:,:,:) = 0._wp 
    254          ! 
    255          DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==! 
     259         zdit(:,:,:) = 0._wp 
     260         zdjt(:,:,:) = 0._wp 
     261         ! 
     262         DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==! 
    256263            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    257264            zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
    258265         END_3D 
    259266         IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level 
    260             DO_2D( 1, 0, 1, 0 )                    ! bottom level 
     267            DO_2D( iij, iij-1, iij, iij-1 )                    ! bottom level 
    261268               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
    262269               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    263270            END_2D 
    264271            IF( ln_isfcav ) THEN                   ! top level (ocean cavities only) 
    265                DO_2D( 1, 0, 1, 0 ) 
     272               DO_2D( iij, iij-1, iij, iij-1 ) 
    266273                  IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 
    267274                  IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 
    … …  
    276283         DO jk = 1, jpkm1 
    277284            !                    !==  Vertical tracer gradient at level jk and jk+1 
    278             DO_2D( 1, 1, 1, 1 ) 
     285            DO_2D( iij, iij, iij, iij ) 
    279286               zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 
    280287            END_2D 
    … …  
    283290            IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1) 
    284291            ELSE 
    285                DO_2D( 1, 1, 1, 1 ) 
     292               DO_2D( iij, iij, iij, iij ) 
    286293                  zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    287294               END_2D 
    … …  
    289296            ! 
    290297            zaei_slp = 0._wp 
     298            zaei_slp_ip1 = 0._wp 
     299            zaei_slp_jp1 = 0._wp 
     300            zaei_slp1 = 0._wp 
    291301            ! 
    292302            IF( ln_botmix_triad ) THEN 
    293                DO ip = 0, 1              !==  Horizontal & vertical fluxes 
    294                   DO kp = 0, 1 
    295                      DO_2D( 1, 0, 1, 0 ) 
    296                         ze1ur = r1_e1u(ji,jj) 
    297                         zdxt  = zdit(ji,jj,jk) * ze1ur 
    298                         ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm) 
    299                         zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr 
    300                         zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    301                         zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp) 
    302                         ! 
    303                         zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 
    304                         ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked.... 
    305                         zah = pahu(ji,jj,jk) 
    306                         zah_slp  = zah * zslope_iso 
    307                         IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew 
    308                         zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
    309                         ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt                 * zbu * ze3wr 
    310                      END_2D 
    311                   END DO 
     303               DO kp = 0, 1              !==  Horizontal & vertical fluxes 
     304                  DO_2D( iij, iij-1, iij, iij-1 ) 
     305                     ze1ur = r1_e1u(ji,jj) 
     306                     zdxt  = zdit(ji,jj,jk) * ze1ur 
     307                     zdxt_ip1  = zdit(ji+1,jj,jk) * r1_e1u(ji+1,jj) 
     308                     ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 
     309                     ze3wr_ip1 = 1._wp / e3w(ji+1,jj,jk+kp,Kmm) 
     310                     zdzt  = zdkt3d(ji,jj,kp) * ze3wr 
     311                     zdzt_ip1  = zdkt3d(ji+1,jj,kp) * ze3wr_ip1 
     312                     ! 
     313                     zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 
     314                     zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm) 
     315                     ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked.... 
     316                     zah = pahu(ji,jj,jk) 
     317                     zah_ip1 = pahu(ji+1,jj,jk) 
     318                     zah_slp  = zah * triadi(ji,jj,jk,1,kp) 
     319                     zah_slp_ip1  = zah_ip1 * triadi(ji+1,jj,jk,1,kp) 
     320                     zah_slp1  = zah * triadi(ji+1,jj,jk,0,kp) 
     321                     IF( ln_ldfeiv )   THEN 
     322                        zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp) 
     323                        zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp) 
     324                        zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp) 
     325                     ENDIF 
     326                     ! round brackets added to fix the order of floating point operations 
     327                     ! needed to ensure halo 1 - halo 2 compatibility 
     328                     zftu(ji   ,jj,jk  ) =  zftu(ji   ,jj,jk )                                                               & 
     329                                         &    - ( ( zah * zdxt + ( zah_slp - zaei_slp ) * zdzt ) * zbu * ze1ur               & 
     330                                         &      + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur  & 
     331                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility 
     332                     ztfw(ji+1,jj,jk+kp) =  ztfw(ji+1,jj,jk+kp)                                                              & 
     333                                         &    - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1              & 
     334                                         &      + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1                           & 
     335                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility 
     336                  END_2D 
    312337               END DO 
    313338               ! 
    314                DO jp = 0, 1 
    315                   DO kp = 0, 1 
    316                      DO_2D( 1, 0, 1, 0 ) 
    317                         ze2vr = r1_e2v(ji,jj) 
    318                         zdyt  = zdjt(ji,jj,jk) * ze2vr 
    319                         ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm) 
    320                         zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr 
    321                         zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    322                         zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
    323                         zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 
    324                         ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked... 
    325                         zah = pahv(ji,jj,jk) 
    326                         zah_slp = zah * zslope_iso 
    327                         IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew 
    328                         zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
    329                         ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt                * zbv * ze3wr 
    330                      END_2D 
    331                   END DO 
     339               DO kp = 0, 1 
     340                  DO_2D( iij, iij-1, iij, iij-1 ) 
     341                     ze2vr = r1_e2v(ji,jj) 
     342                     zdyt  = zdjt(ji,jj,jk) * ze2vr 
     343                     zdyt_jp1  = zdjt(ji,jj+1,jk) * r1_e2v(ji,jj+1) 
     344                     ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 
     345                     ze3wr_jp1 = 1._wp / e3w(ji,jj+1,jk+kp,Kmm) 
     346                     zdzt  = zdkt3d(ji,jj,kp) * ze3wr 
     347                     zdzt_jp1  = zdkt3d(ji,jj+1,kp) * ze3wr_jp1 
     348                     zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     349                     zbv_jp1 = 0.25_wp * e1e2v(ji,jj+1) * e3v(ji,jj+1,jk,Kmm) 
     350                     ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked.... 
     351                     zah = pahv(ji,jj,jk)          ! pahv(ji,jj+jp,jk)  ???? 
     352                     zah_jp1 = pahv(ji,jj+1,jk) 
     353                     zah_slp = zah * triadj(ji,jj,jk,1,kp) 
     354                     zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp) 
     355                     zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp) 
     356                     IF( ln_ldfeiv )   THEN 
     357                        zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 
     358                        zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp) 
     359                        zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp) 
     360                     ENDIF 
     361                     ! round brackets added to fix the order of floating point operations 
     362                     ! needed to ensure halo 1 - halo 2 compatibility 
     363                     zftv(ji,jj  ,jk   ) =  zftv(ji,jj  ,jk   )                                                              & 
     364                                         &    - ( ( zah * zdyt + ( zah_slp - zaei_slp ) * zdzt ) * zbv * ze2vr               & 
     365                                         &      + ( zah * zdyt + zah_slp1 * zdzt_jp1 - zaei_slp1 * zdzt_jp1 ) * zbv * ze2vr  & 
     366                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility 
     367                     ztfw(ji,jj+1,jk+kp) =  ztfw(ji,jj+1,jk+kp)                                                              & 
     368                                         &    - ( ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1             & 
     369                                         &      + ( zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1                           & 
     370                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility 
     371                  END_2D 
    332372               END DO 
    333373               ! 
    334374            ELSE 
    335375               ! 
    336                DO ip = 0, 1               !==  Horizontal & vertical fluxes 
    337                   DO kp = 0, 1 
    338                      DO_2D( 1, 0, 1, 0 ) 
    339                         ze1ur = r1_e1u(ji,jj) 
    340                         zdxt  = zdit(ji,jj,jk) * ze1ur 
    341                         ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm) 
    342                         zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr 
    343                         zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 
    344                         zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp) 
    345                         ! 
    346                         zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 
    347                         ! ln_botmix_triad is .F. mask zah for bottom half cells 
    348                         zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ???? 
    349                         zah_slp  = zah * zslope_iso 
    350                         IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! aeit(ji+ip,jj,jk)*zslope_skew 
    351                         zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
    352                         ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 
    353                      END_2D 
    354                   END DO 
     376               DO kp = 0, 1               !==  Horizontal & vertical fluxes 
     377                  DO_2D( iij, iij-1, iij, iij-1 ) 
     378                     ze1ur = r1_e1u(ji,jj) 
     379                     zdxt  = zdit(ji,jj,jk) * ze1ur 
     380                     zdxt_ip1  = zdit(ji+1,jj,jk) * r1_e1u(ji+1,jj) 
     381                     ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 
     382                     ze3wr_ip1 = 1._wp / e3w(ji+1,jj,jk+kp,Kmm) 
     383                     zdzt  = zdkt3d(ji,jj,kp) * ze3wr 
     384                     zdzt_ip1  = zdkt3d(ji+1,jj,kp) * ze3wr_ip1 
     385                     ! 
     386                     zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 
     387                     zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm) 
     388                     ! ln_botmix_triad is .F. mask zah for bottom half cells 
     389                     zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ???? 
     390                     zah_ip1 = pahu(ji+1,jj,jk) * umask(ji+1,jj,jk+kp) 
     391                     zah_slp  = zah * triadi(ji,jj,jk,1,kp) 
     392                     zah_slp_ip1  = zah_ip1 * triadi(ji+1,jj,jk,1,kp) 
     393                     zah_slp1  = zah * triadi(ji+1,jj,jk,0,kp) 
     394                     IF( ln_ldfeiv )   THEN 
     395                        zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp) 
     396                        zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp) 
     397                        zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp) 
     398                     ENDIF 
     399                     ! round brackets added to fix the order of floating point operations 
     400                     ! needed to ensure halo 1 - halo 2 compatibility 
     401                     zftu(ji   ,jj,jk  ) =  zftu(ji   ,jj,jk )                                                               & 
     402                                         &    - ( ( zah * zdxt + ( zah_slp - zaei_slp ) * zdzt ) * zbu * ze1ur               & 
     403                                         &      + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur  & 
     404                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility 
     405                     ztfw(ji+1,jj,jk+kp) =  ztfw(ji+1,jj,jk+kp)                                                              & 
     406                                         &    - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1              & 
     407                                         &      + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1                           & 
     408                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility 
     409                  END_2D 
    355410               END DO 
    356411               ! 
    357                DO jp = 0, 1 
    358                   DO kp = 0, 1 
    359                      DO_2D( 1, 0, 1, 0 ) 
    360                         ze2vr = r1_e2v(ji,jj) 
    361                         zdyt  = zdjt(ji,jj,jk) * ze2vr 
    362                         ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm) 
    363                         zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr 
    364                         zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 
    365                         zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp) 
    366                         zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 
    367                         ! ln_botmix_triad is .F. mask zah for bottom half cells 
    368                         zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ???? 
    369                         zah_slp = zah * zslope_iso 
    370                         IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! aeit(ji,jj+jp,jk)*zslope_skew 
    371                         zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
    372                         ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 
    373                      END_2D 
    374                   END DO 
     412               DO kp = 0, 1 
     413                  DO_2D( iij, iij-1, iij, iij-1 ) 
     414                     ze2vr = r1_e2v(ji,jj) 
     415                     zdyt  = zdjt(ji,jj,jk) * ze2vr 
     416                     zdyt_jp1  = zdjt(ji,jj+1,jk) * r1_e2v(ji,jj+1) 
     417                     ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 
     418                     ze3wr_jp1 = 1._wp / e3w(ji,jj+1,jk+kp,Kmm) 
     419                     zdzt  = zdkt3d(ji,jj,kp) * ze3wr 
     420                     zdzt_jp1  = zdkt3d(ji,jj+1,kp) * ze3wr_jp1 
     421                     zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 
     422                     zbv_jp1 = 0.25_wp * e1e2v(ji,jj+1) * e3v(ji,jj+1,jk,Kmm) 
     423                     ! ln_botmix_triad is .F. mask zah for bottom half cells 
     424                     zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ???? 
     425                     zah_jp1 = pahv(ji,jj+1,jk) * vmask(ji,jj+1,jk+kp) 
     426                     zah_slp = zah * triadj(ji,jj,jk,1,kp) 
     427                     zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp) 
     428                     zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp) 
     429                     IF( ln_ldfeiv )   THEN 
     430                        zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 
     431                        zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp) 
     432                        zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp) 
     433                     ENDIF 
     434                     ! round brackets added to fix the order of floating point operations 
     435                     ! needed to ensure halo 1 - halo 2 compatibility 
     436                     zftv(ji,jj  ,jk   ) =  zftv(ji,jj  ,jk   )                                                              & 
     437                                         &    - ( ( zah * zdyt + ( zah_slp - zaei_slp ) * zdzt ) * zbv * ze2vr               & 
     438                                         &      + ( zah * zdyt + zah_slp1 * zdzt_jp1 - zaei_slp1 * zdzt_jp1 ) * zbv * ze2vr  & 
     439                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility 
     440                     ztfw(ji,jj+1,jk+kp) =  ztfw(ji,jj+1,jk+kp)                                                              & 
     441                                         &    - ( ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1             & 
     442                                         &      + ( zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1                           & 
     443                                         &      )                                                                            ! bracket for halo 1 - halo 2 compatibility 
     444                  END_2D 
    375445               END DO 
    376446            ENDIF 
    377447            !                             !==  horizontal divergence and add to the general trend  ==! 
    378             DO_2D( 0, 0, 0, 0 ) 
    379                pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    & 
    380                   &                       + zsign * (  zftu(ji-1,jj  ,jk) - zftu(ji,jj,jk)       & 
    381                   &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   & 
    382                   &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  ) 
     448            DO_2D( iij-1, iij-1, iij-1, iij-1 ) 
     449               ! round brackets added to fix the order of floating point operations 
     450               ! needed to ensure halo 1 - halo 2 compatibility 
     451               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)                                                & 
     452                  &                       + zsign * ( ( zftu(ji-1,jj  ,jk) - zftu(ji,jj,jk)             & 
     453                  &                                   )                                                 & ! bracket for halo 1 - halo 2 compatibility 
     454                  &                                 + ( zftv(ji,jj-1,jk) - zftv(ji,jj,jk)               & 
     455                  &                                   )                                                 & ! bracket for halo 1 - halo 2 compatibility 
     456                  &                                 ) / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  ) 
    383457            END_2D 
    384458            ! 
    … …  
    387461         !                                !==  add the vertical 33 flux  ==! 
    388462         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
    389             DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 
     463            DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 
    390464               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   & 
    391465                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    … …  
    395469            SELECT CASE( kpass ) 
    396470            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    397                DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 
     471               DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 
    398472                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             & 
    399473                     &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    400474               END_3D 
    401475            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp. 
    402                DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 
     476               DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    403477                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      & 
    404478                     &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   & 
    … …  
    408482         ENDIF 
    409483         ! 
    410          DO_3D( 0, 0, 0, 0, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==! 
     484         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==! 
    411485            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    & 
    412486            &                                  + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/tramle.F90

    r14644 r14986  
    8787      INTEGER                     , INTENT(in   ) ::   Kmm        ! ocean time level index 
    8888      CHARACTER(len=3)            , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    89       REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components 
    90       REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   pv         ! out: same 3  transport components 
    91       REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport 
     89      ! TEMP: [tiling] Can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 
     90      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components 
     91      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components 
     92      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport 
    9293      ! 
    9394      INTEGER  ::   ji, jj, jk          ! dummy loop indices 
    … …  
    9697      REAL(wp) ::   zcvw, zmvw          !   -      - 
    9798      INTEGER , DIMENSION(A2D(nn_hls))     :: inml_mle 
    98       REAL(wp), DIMENSION(A2D(nn_hls))     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_MH 
     99      REAL(wp), DIMENSION(A2D(nn_hls))     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 
    99100      REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpsi_uw, zpsi_vw 
    100       ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support) 
    101       REAL(wp), DIMENSION(:,:),   ALLOCATABLE, SAVE :: zLf_NH 
    102       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zpsiu_mle, zpsiv_mle 
    103101      !!---------------------------------------------------------------------- 
    104102      ! 
    … …  
    110108         SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
    111109         CASE ( 0 )                                               != min of the 2 neighbour MLDs 
    112             DO_2D( 1, 0, 1, 0 ) 
     110            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    113111               zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 
    114112               zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 
    115113            END_2D 
    116114         CASE ( 1 )                                               != average of the 2 neighbour MLDs 
    117             DO_2D( 1, 0, 1, 0 ) 
     115            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    118116               zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 
    119117               zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 
    120118            END_2D 
    121119         CASE ( 2 )                                               != max of the 2 neighbour MLDs 
    122             DO_2D( 1, 0, 1, 0 ) 
     120            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    123121               zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 
    124122               zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 
    … …  
    126124         END SELECT 
    127125         IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
    128             DO_2D( 1, 0, 1, 0 ) 
     126            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    129127               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            & 
    130128                    &           * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
    … …  
    137135            ! 
    138136         ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
    139             DO_2D( 1, 0, 1, 0 ) 
     137            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    140138               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj)               & 
    141139                    &                  * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) 
    … …  
    149147         !                                      !==  MLD used for MLE  ==! 
    150148         !                                                ! compute from the 10m density to deal with the diurnal cycle 
    151          DO_2D( 1, 1, 1, 1 ) 
     149         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    152150            inml_mle(ji,jj) = mbkt(ji,jj) + 1                    ! init. to number of ocean w-level (T-level + 1) 
    153151         END_2D 
    154152         IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m 
    155            DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m) 
     153           DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m) 
    156154              IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer 
    157155           END_3D 
    … …  
    163161         zbm (:,:) = 0._wp 
    164162         zn2 (:,:) = 0._wp 
    165          DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer 
     163         DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer 
    166164            zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points 
    167165            zmld(ji,jj) = zmld(ji,jj) + zc 
    … …  
    172170         SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts 
    173171         CASE ( 0 )                                               != min of the 2 neighbour MLDs 
    174             DO_2D( 1, 0, 1, 0 ) 
     172            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    175173               zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 
    176174               zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 
    177175            END_2D 
    178176         CASE ( 1 )                                               != average of the 2 neighbour MLDs 
    179             DO_2D( 1, 0, 1, 0 ) 
     177            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    180178               zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 
    181179               zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 
    182180            END_2D 
    183181         CASE ( 2 )                                               != max of the 2 neighbour MLDs 
    184             DO_2D( 1, 0, 1, 0 ) 
     182            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    185183               zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 
    186184               zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) 
    … …  
    188186         END SELECT 
    189187         !                                                ! convert density into buoyancy 
    190          DO_2D( 1, 1, 1, 1 ) 
     188         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    191189            zbm(ji,jj) = + grav * zbm(ji,jj) / MAX( e3t(ji,jj,1,Kmm), zmld(ji,jj) ) 
    192190         END_2D 
    … …  
    201199         ! 
    202200         IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation 
    203             DO_2D( 1, 0, 1, 0 ) 
     201            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    204202               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            & 
    205203                    &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   & 
    … …  
    212210            ! 
    213211         ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 
    214             DO_2D( 1, 0, 1, 0 ) 
     212            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    215213               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               & 
    216214                    &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) 
    … …  
    222220         ! 
    223221         IF( nn_conv == 1 ) THEN              ! No MLE in case of convection 
    224             DO_2D( 1, 0, 1, 0 ) 
     222            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    225223               IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp 
    226224               IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp 
    … …  
    230228      ENDIF  ! end of ln_osm_mle conditional 
    231229    !                                      !==  structure function value at uw- and vw-points  ==! 
    232     DO_2D( 1, 0, 1, 0 ) 
     230    DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    233231       zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall)                   ! hu --> 1/hu 
    234232       zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall)  
    … …  
    238236    zpsi_vw(:,:,:) = 0._wp 
    239237    ! 
    240       DO_3D( 1, 0, 1, 0, 2, ikmax )                ! start from 2 : surface value = 0 
     238      DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, ikmax )                ! start from 2 : surface value = 0 
     239       
    241240         zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj) 
    242241         zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj) 
    … …  
    252251      !                                      !==  transport increased by the MLE induced transport ==! 
    253252      DO jk = 1, ikmax 
    254          DO_2D( 1, 0, 1, 0 )                      ! CAUTION pu,pv must be defined at row/column i=1 / j=1 
     253         DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    255254            pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 
    256255            pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 
    257256         END_2D 
    258          DO_2D( 0, 0, 0, 0 ) 
     257         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
    259258            pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   & 
    260259               &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) * wmask(ji,jj,1) 
    … …  
    262261      END DO 
    263262 
    264       ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support) 
    265263      IF( cdtype == 'TRA') THEN              !==  outputs  ==! 
    266          IF( ntile == 0 .OR. ntile == 1 ) THEN                             ! Do only on the first tile 
    267             ALLOCATE( zLf_NH(jpi,jpj), zpsiu_mle(jpi,jpj,jpk), zpsiv_mle(jpi,jpj,jpk) ) 
    268             zpsiu_mle(:,:,:) = 0._wp ; zpsiv_mle(:,:,:) = 0._wp 
    269          ENDIF 
    270264         ! 
    271265         IF (ln_osm_mle.and.ln_zdfosm) THEN 
    … …  
    279273         ENDIF 
    280274         ! 
     275         CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f 
     276         ! 
    281277         ! divide by cross distance to give streamfunction with dimensions m^2/s 
    282278         DO_3D( 0, 0, 0, 0, 1, ikmax+1 ) 
    283             zpsiu_mle(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj) 
    284             zpsiv_mle(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj) 
     279            zpsi_uw(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj) 
     280            zpsi_vw(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj) 
    285281         END_3D 
    286  
    287          IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile 
    288             CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f 
    289             CALL iom_put( "psiu_mle", zpsiu_mle )    ! i-mle streamfunction 
    290             CALL iom_put( "psiv_mle", zpsiv_mle )    ! j-mle streamfunction 
    291             DEALLOCATE( zLf_NH, zpsiu_mle, zpsiv_mle ) 
    292          ENDIF 
     282         CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction 
     283         CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction 
    293284      ENDIF 
    294285      ! 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/tranpc.F90

    r14644 r14986  
    1717   USE oce            ! ocean dynamics and active tracers 
    1818   USE dom_oce        ! ocean space and time domain 
    19    ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed) 
    20    USE domtile 
    2119   USE phycst         ! physical constants 
    2220   USE zdf_oce        ! ocean vertical physics 
    … …  
    8280      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is 
    8381      INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1" 
    84       INTEGER :: isi, isj, iei, iej 
    8582      LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu" 
    8683      !!---------------------------------------------------------------------- 
    … …  
    106103         CALL bn2    ( CASTWP(pts(:,:,:,:,Kaa)), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points) 
    107104         ! 
    108          IF( ntile == 0 .OR. ntile == 1 ) nnpcc = 0         ! Do only on the first tile 
    109          ! 
    110          IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling 
    111          IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF 
    112          IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF 
    113          IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF 
    114          ! 
    115          DO_2D( isi, iei, isj, iej )                        ! interior column only 
     105         IF( .NOT. l_istiled .OR. ntile == 1 ) nnpcc = 0         ! Do only on the first tile 
     106         ! 
     107         DO_2D_OVR( 0, 0, 0, 0 )                        ! interior column only 
    116108            ! 
    117109            IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points 
    … …  
    320312         ENDIF 
    321313         ! 
    322          IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain 
     314         IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only for the full domain 
    323315            IF( lwp .AND. l_LB_debug ) THEN 
    324316               WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', nnpcc 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/traqsr.F90

    r14644 r14986  
    109109      ! 
    110110      INTEGER  ::   ji, jj, jk               ! dummy loop indices 
    111       INTEGER  ::   irgb, isi, iei, isj, iej ! local integers 
     111      INTEGER  ::   irgb                    ! local integers 
    112112      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars 
    113113      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         - 
    … …  
    122122      IF( ln_timing )   CALL timing_start('tra_qsr') 
    123123      ! 
    124       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     124      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    125125         IF( kt == nit000 ) THEN 
    126126            IF(lwp) WRITE(numout,*) 
    … …  
    138138      !                         !  before qsr induced heat content  ! 
    139139      !                         !-----------------------------------! 
    140       IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling 
    141       IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF 
    142       IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF 
    143       IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF 
    144  
    145140      IF( kt == nit000 ) THEN          !==  1st time step  ==! 
    146141         IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN    ! read in restart 
    147142            z1_2 = 0.5_wp 
    148             IF( ntile == 0 .OR. ntile == 1 )  THEN                        ! Do only on the first tile 
     143            IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                        ! Do only on the first tile 
    149144               IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file' 
    150145               CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b )   ! before heat content trend due to Qsr flux 
    … …  
    152147         ELSE                                           ! No restart or Euler forward at 1st time step 
    153148            z1_2 = 1._wp 
    154             DO_3D( isi, iei, isj, iej, 1, jpk ) 
     149            DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 
    155150               qsr_hc_b(ji,jj,jk) = 0._wp 
    156151            END_3D 
    … …  
    158153      ELSE                             !==  Swap of qsr heat content  ==! 
    159154         z1_2 = 0.5_wp 
    160          DO_3D( isi, iei, isj, iej, 1, jpk ) 
     155         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 
    161156            qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk) 
    162157         END_3D 
    … …  
    169164      CASE( np_BIO )                   !==  bio-model fluxes  ==! 
    170165         ! 
    171          DO_3D( isi, iei, isj, iej, 1, nksr ) 
     166         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr ) 
    172167            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) 
    173168         END_3D 
    … …  
    180175         ! 
    181176         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll 
    182             IF( ntile == 0 .OR. ntile == 1 )  THEN                                         ! Do only for the full domain 
    183                IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )            ! Use full domain 
     177            IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                                         ! Do only for the full domain 
     178               IF( ln_tile ) CALL dom_tile_stop( ldhold=.TRUE. )             ! Use full domain 
    184179               CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step 
    185                IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 )            ! Revert to tile domain 
     180               IF( ln_tile ) CALL dom_tile_start( ldhold=.TRUE. )            ! Revert to tile domain 
    186181            ENDIF 
    187182            ! 
    … …  
    191186            ! most expensive calculations) 
    192187            ! 
    193             DO_2D( isi, iei, isj, iej ) 
     188            DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
    194189                       ! zlogc = log(zchl) 
    195190               zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) ) 
    … …  
    210205 
    211206! 
    212             DO_3D( isi, iei, isj, iej, 1, nksr + 1 ) 
     207            DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr + 1 ) 
    213208               ! zchl    = ALOG( ze0(ji,jj) ) 
    214209               zlogc = ze0(ji,jj) 
    … …  
    240235         ! 
    241236         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B 
    242          DO_2D( isi, iei, isj, iej ) 
     237         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
    243238            ze0(ji,jj) = rn_abs * qsr(ji,jj) 
    244239            ze1(ji,jj) = zcoef  * qsr(ji,jj) 
    … …  
    251246         ! 
    252247         !                                    !* interior equi-partition in R-G-B depending on vertical profile of Chl 
    253          DO_3D( isi, iei, isj, iej, 2, nksr + 1 ) 
     248         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr + 1 ) 
    254249            ze3t = e3t(ji,jj,jk-1,Kmm) 
    255250            irgb = NINT( ztmp3d(ji,jj,jk) ) 
    … …  
    265260         END_3D 
    266261         ! 
    267          DO_3D( isi, iei, isj, iej, 1, nksr )          !* now qsr induced heat content 
     262         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr )          !* now qsr induced heat content 
    268263            qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 
    269264         END_3D 
    … …  
    275270         zz0 =        rn_abs   * r1_rho0_rcp      ! surface equi-partition in 2-bands 
    276271         zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 
    277          DO_3D( isi, iei, isj, iej, 1, nksr )          !* now qsr induced heat content 
     272         DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr )          !* now qsr induced heat content 
    278273            zc0 = zz0 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk  ,Kmm)*xsi1r ) 
    279274            zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r ) 
    … …  
    293288      ! 
    294289      ! sea-ice: store the 1st ocean level attenuation coefficient 
    295       DO_2D( isi, iei, isj, iej ) 
     290      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
    296291         IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 
    297292         ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp 
    … …  
    299294      END_2D 
    300295      ! 
    301       ! TEMP: [tiling] This change not necessary and working array can use A2D(nn_hls) if using XIOS (subdomain support) 
    302       IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain 
    303          IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution 
    304             ALLOCATE( zetot(jpi,jpj,jpk) ) 
    305             zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero 
    306             DO jk = nksr, 1, -1 
    307                zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 
    308             END DO 
    309             CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation 
    310             DEALLOCATE( zetot ) 
    311          ENDIF 
    312       ENDIF 
    313       ! 
    314       IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile 
     296      IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution 
     297         ALLOCATE( zetot(A2D(nn_hls),jpk) ) 
     298         zetot(:,:,nksr+1:jpk) = 0._wp     ! below ~400m set to zero 
     299         DO_3DS(0, 0, 0, 0, nksr, 1, -1) 
     300            zetot(ji,jj,jk) = zetot(ji,jj,jk+1) + qsr_hc(ji,jj,jk) * rho0_rcp 
     301         END_3D 
     302         CALL iom_put( 'qsr3d', zetot )   ! 3D distribution of shortwave Radiation 
     303         DEALLOCATE( zetot ) 
     304      ENDIF 
     305      ! 
     306      IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile 
    315307         IF( lrst_oce ) THEN     ! write in the ocean restart file 
    316308            CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      ) 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/trasbc.F90

    r14644 r14986  
    7878      ! 
    7979      INTEGER  ::   ji, jj, jk, jn               ! dummy loop indices 
    80       INTEGER  ::   ikt, ikb, isi, iei, isj, iej ! local integers 
     80      INTEGER  ::   ikt, ikb                    ! local integers 
    8181      REAL(wp) ::   zfact, z1_e3t, zdep, ztim    ! local scalar 
    8282      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    … …  
    8585      IF( ln_timing )   CALL timing_start('tra_sbc') 
    8686      ! 
    87       IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     87      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    8888         IF( kt == nit000 ) THEN 
    8989            IF(lwp) WRITE(numout,*) 
    … …  
    9999      ENDIF 
    100100      ! 
    101       IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling 
    102       IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF 
    103       IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF 
    104       IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF 
    105  
    106101!!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist) 
    107102      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration 
    108          DO_2D( isi, iei, isj, iej ) 
     103         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
    109104            qns(ji,jj) = qns(ji,jj) + qsr(ji,jj)      ! total heat flux in qns 
    110105            qsr(ji,jj) = 0._wp                        ! qsr set to zero 
    … …  
    119114         IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN      ! Restart: read in restart file 
    120115            zfact = 0.5_wp 
    121             IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
     116            IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile 
    122117               IF(lwp) WRITE(numout,*) '          nit000-1 sbc tracer content field read in the restart file' 
    123118               sbc_tsc(:,:,:) = 0._wp 
    … …  
    127122         ELSE                                             ! No restart or restart not found: Euler forward time stepping 
    128123            zfact = 1._wp 
    129             DO_2D( isi, iei, isj, iej ) 
     124            DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
    130125               sbc_tsc(ji,jj,:) = 0._wp 
    131126               sbc_tsc_b(ji,jj,:) = 0._wp 
    … …  
    134129      ELSE                                !* other time-steps: swap of forcing fields 
    135130         zfact = 0.5_wp 
    136          DO_2D( isi, iei, isj, iej ) 
     131         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
    137132            sbc_tsc_b(ji,jj,:) = sbc_tsc(ji,jj,:) 
    138133         END_2D 
    139134      ENDIF 
    140135      !                             !==  Now sbc tracer content fields  ==! 
    141       DO_2D( isi, iei, isj, iej ) 
     136      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 
    142137         sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj)   ! non solar heat flux 
    143138         sbc_tsc(ji,jj,jp_sal) = r1_rho0     * sfx(ji,jj)   ! salt flux due to freezing/melting 
    144139      END_2D 
    145140      IF( ln_linssh ) THEN                !* linear free surface 
    146          DO_2D( isi, iei, isj, iej )                    !==>> add concentration/dilution effect due to constant volume cell 
     141         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )                    !==>> add concentration/dilution effect due to constant volume cell 
    147142            sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm) 
    148143            sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm) 
    149144         END_2D                                 !==>> output c./d. term 
    150          IF( ntile == 0 .OR. ntile == nijtile )  THEN             ! Do only on the last tile 
    151             IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) 
    152             IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) ) 
    153          ENDIF 
     145         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) 
     146         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) ) 
    154147      ENDIF 
    155148      ! 
    … …  
    161154      END DO 
    162155      ! 
    163       IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile 
     156      IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only on the last tile 
    164157         IF( lrst_oce ) THEN           !==  write sbc_tsc in the ocean restart file  ==! 
    165158            CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) ) 
    … …  
    187180      ENDIF 
    188181 
    189       IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile 
    190          IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) )   ! runoff term on sst 
    191          IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) )   ! runoff term on sss 
    192       ENDIF 
     182      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) )   ! runoff term on sst 
     183      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) )   ! runoff term on sss 
    193184 
    194185#if defined key_asminc 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/trazdf.F90

    r14644 r14986  
    6565      ! 
    6666      IF( kt == nit000 )  THEN 
    67          IF( ntile == 0 .OR. ntile == 1 )  THEN                   ! Do only on the first tile 
     67         IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                   ! Do only on the first tile 
    6868            IF(lwp)WRITE(numout,*) 
    6969            IF(lwp)WRITE(numout,*) 'tra_zdf : implicit vertical mixing on T & S' 

  • NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/TRA/zpshde.F90

    r14930 r14986  
    4747      INTEGER                     , INTENT(in   )           ::  Kmm         ! ocean time level index 
    4848      INTEGER                     , INTENT(in   )           ::  kjpt        ! number of tracers 
    49       REAL(dp), DIMENSION(:,:,:,:), INTENT(inout)           ::  pta         ! 4D tracers fields 
     49      REAL(dp), DIMENSION(:,:,:,:), INTENT(in   )           ::  pta         ! 4D tracers fields 
    5050      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts 
    51       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout), OPTIONAL ::  prd         ! 3D density anomaly fields 
     51      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    5252      REAL(wp), DIMENSION(:,:)    , INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom) 
    5353      ! 
    … …  
    111111      INTEGER                                , INTENT(in   )           ::  kjpt        ! number of tracers 
    112112      INTEGER                                , INTENT(in   )           ::  ktta, ktgt, ktrd, ktgr 
    113       REAL(dp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(inout)           ::  pta         ! 4D tracers fields 
     113      REAL(dp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in   )           ::  pta         ! 4D tracers fields 
    114114      REAL(wp), DIMENSION(A2D_T(ktgt)    ,KJPT), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts 
    115       REAL(wp), DIMENSION(A2D_T(ktrd),JPK     ), INTENT(inout), OPTIONAL ::  prd         ! 3D density anomaly fields 
     115      REAL(wp), DIMENSION(A2D_T(ktrd),JPK     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    116116      REAL(wp), DIMENSION(A2D_T(ktgr)         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom) 
    117117      ! 
    … …  
    124124      ! 
    125125      IF( ln_timing )   CALL timing_start( 'zps_hde') 
    126       IF (nn_hls.EQ.2) THEN 
    127          CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp) 
    128          IF(PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp) 
    129       END IF 
    130126      ! 
    131127      pgtu(:,:,:) = 0._wp   ;   zti (:,:,:) = 0._wp   ;   zhi (:,:) = 0._wp 
    … …  
    134130      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    135131         ! 
    136          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )              ! Gradient of density at the last level 
     132         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )              ! Gradient of density at the last level 
    137133            iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    138134            ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
    … …  
    173169      END DO 
    174170      ! 
    175       IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond. 
     171      IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond. 
    176172      ! 
    177173      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
    … …  
    206202            ENDIF 
    207203         END_2D 
    208          IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions 
     204         IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions 
    209205         ! 
    210206      END IF 
    … …  
    221217      INTEGER                     , INTENT(in   )           ::  Kmm          ! ocean time level index 
    222218      INTEGER                     , INTENT(in   )           ::  kjpt         ! number of tracers 
    223       REAL(dp), DIMENSION(:,:,:,:), INTENT(inout)           ::  pta          ! 4D tracers fields 
     219      REAL(dp), DIMENSION(:,:,:,:), INTENT(in   )           ::  pta          ! 4D tracers fields 
    224220      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts 
    225221      REAL(wp), DIMENSION(:,:,:)  , INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 
    226       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout), OPTIONAL ::  prd          ! 3D density anomaly fields 
     222      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ), OPTIONAL ::  prd          ! 3D density anomaly fields 
    227223      REAL(wp), DIMENSION(:,:)    , INTENT(  out), OPTIONAL ::  pgru, pgrv   ! hor. grad of prd at u- & v-pts (bottom) 
    228224      REAL(wp), DIMENSION(:,:)    , INTENT(  out), OPTIONAL ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 
    … …  
    291287      INTEGER                                , INTENT(in   )           ::  kjpt         ! number of tracers 
    292288      INTEGER                                , INTENT(in   )           ::  ktta, ktgt, ktgti, ktrd, ktgr, ktgri 
    293       REAL(dp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(inout)           ::  pta          ! 4D tracers fields 
     289      REAL(dp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in   )           ::  pta          ! 4D tracers fields 
    294290      REAL(wp), DIMENSION(A2D_T(ktgt)    ,KJPT), INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts 
    295291      REAL(wp), DIMENSION(A2D_T(ktgti)   ,KJPT), INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 
    296       REAL(wp), DIMENSION(A2D_T(ktrd),JPK     ), INTENT(inout), OPTIONAL ::  prd          ! 3D density anomaly fields 
     292      REAL(wp), DIMENSION(A2D_T(ktrd),JPK     ), INTENT(in   ), OPTIONAL ::  prd          ! 3D density anomaly fields 
    297293      REAL(wp), DIMENSION(A2D_T(ktgr)         ), INTENT(  out), OPTIONAL ::  pgru, pgrv   ! hor. grad of prd at u- & v-pts (bottom) 
    298294      REAL(wp), DIMENSION(A2D_T(ktgri)        ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 
    … …  
    307303      IF( ln_timing )   CALL timing_start( 'zps_hde_isf') 
    308304      ! 
    309       IF (nn_hls.EQ.2) THEN 
    310          CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp) 
    311          IF (PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp) 
    312       END IF 
    313  
    314305      pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp 
    315306      pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp 
    … …  
    319310      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    320311         ! 
    321          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     312         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    322313 
    323314            iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    … …  
    359350      END DO 
    360351      ! 
    361       IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond. 
     352      IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond. 
    362353 
    363354      ! horizontal derivative of density anomalies (rd) 
    … …  
    401392         END_2D 
    402393 
    403          IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions 
     394         IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp )   ! Lateral boundary conditions 
    404395         ! 
    405396      END IF 
    … …  
    408399      ! 
    409400      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
    410          DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     401         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    411402            iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1 
    412403            ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1 
    … …  
    452443         ! 
    453444      END DO 
    454       IF (nn_hls.EQ.1) CALL lbc_lnk( 'zpshde', pgtui(:,:,:), 'U', -1.0_wp , pgtvi(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond. 
     445      IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgtui(:,:,:), 'U', -1.0_wp , pgtvi(:,:,:), 'V', -1.0_wp )   ! Lateral boundary cond. 
    455446 
    456447      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
    … …  
    491482 
    492483         END_2D 
    493          IF (nn_hls.EQ.1) THEN 
    494            CALL lbc_lnk( 'zpshde', pgrui, 'U', -1.0_wp) 
     484         IF (nn_hls==1) THEN 
     485           CALL lbc_lnk( 'zpshde', pgrui, 'U', -1.0_wp )   ! Lateral boundary conditions 
    495486           CALL lbc_lnk( 'zpshde', pgrvi, 'V', -1.0_wp )   ! Lateral boundary conditions 
    496487         ENDIF 
Note: See TracChangeset for help on using the changeset viewer.