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 14631 for NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_triad.F90 – NEMO

Ignore:
Timestamp:
2021-03-23T17:55:26+01:00 (3 years ago)
Author:
hadcv
Message:

#2600: Update TRA tiling after [14576] merge (includes tra_ldf_triad changes missing from [14607])

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_triad.F90

    r14576 r14631  
    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 
     
    111109      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
    112110      INTEGER  ::  ip,jp,kp         ! dummy loop indices 
    113       INTEGER  ::  ierr            ! local integer 
     111      ! NOTE: [halo1-halo2] 
     112      INTEGER  ::  ierr, iij        ! local integer 
    114113      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3    ! local scalars 
    115114      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4    !   -      - 
     
    119118      REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt 
    120119      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 
     120      REAL(wp), DIMENSION(A2D(nn_hls),0:1) ::   zdkt3d                                           ! vertical tracer gradient at 2 levels 
     121      REAL(wp), DIMENSION(A2D(nn_hls)    ) ::   z2d                                              ! 2D workspace 
     122      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zdit, zdjt, zpsi_uw, zpsi_vw   ! 3D     - 
     123      ! NOTE: [halo1-halo2] 
     124      REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: zftu, zftv 
     125      REAL(wp), DIMENSION(A2D(nn_hls),jpk,2,2) :: ztfw 
    126126      !!---------------------------------------------------------------------- 
    127127      ! 
     
    142142      ENDIF 
    143143      ! 
     144      ! NOTE: [halo1-halo2] 
     145      ! Define pt_rhs halo points for multi-point haloes in bilaplacian case 
     146      IF( nldf_tra == np_blp_it .AND. kpass == 1 ) THEN ; iij = nn_hls 
     147      ELSE                                              ; iij = 1 
     148      ENDIF 
     149 
     150      ! 
    144151      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0) 
    145152      ELSE                    ;   zsign = -1._wp 
     
    153160         ! 
    154161         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpk ) 
    155          DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )  
     162         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
    156163            akz     (ji,jj,jk) = 0._wp 
    157164            ah_wslp2(ji,jj,jk) = 0._wp 
    158165         END_3D 
     166         ! 
     167         ! NOTE: [halo1-halo2] 
     168         zftu(:,:,:,:) = 0._wp 
     169         zftv(:,:,:,:) = 0._wp 
    159170         ! 
    160171         DO ip = 0, 1                            ! i-k triads 
     
    169180                  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) 
    170181                  zslope2 = zslope2 *zslope2 
    171                   ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 
    172                   akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + zah * r1_e1u(ji-ip,jj)       & 
    173                      &                                                      * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 
     182                  ! NOTE: [halo1-halo2] 
     183                  zftu(ji,jj,jk+kp,ip+1) = zftu(ji,jj,jk+kp,ip+1) + & 
     184                     &                     zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 
     185                  zftv(ji,jj,jk+kp,ip+1) = zftv(ji,jj,jk+kp,ip+1) + & 
     186                     &                     zah * r1_e1u(ji-ip,jj) * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 
    174187                     ! 
    175188               END_3D 
    176189            END DO 
    177190         END DO 
     191         ! 
     192         ! NOTE: [halo1-halo2] Use DO_3D instead of DO_3D_OVR 
     193         ! ip loop contributions added here to ensure consistent floating point arithmetic for different nn_hls 
     194         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
     195            ah_wslp2(ji,jj,jk) = ah_wslp2(ji,jj,jk) + (zftu(ji,jj,jk,1) + zftu(ji,jj,jk,2)) 
     196            akz     (ji,jj,jk) = akz     (ji,jj,jk) + (zftv(ji,jj,jk,1) + zftv(ji,jj,jk,2)) 
     197         END_3D 
     198         ! 
     199         ! NOTE: [halo1-halo2] 
     200         zftu(:,:,:,:) = 0._wp 
     201         zftv(:,:,:,:) = 0._wp 
    178202         ! 
    179203         DO jp = 0, 1                            ! j-k triads 
     
    189213                  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) 
    190214                  zslope2 = zslope2 * zslope2 
    191                   ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 
    192                   akz     (ji,jj,jk+kp) = akz     (ji,jj,jk+kp) + zah * r1_e2v(ji,jj-jp)     & 
    193                      &                                                      * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 
     215                  ! NOTE: [halo1-halo2] 
     216                  zftu(ji,jj,jk+kp,jp+1) = zftu(ji,jj,jk+kp,jp+1) + & 
     217                     &                     zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 
     218                  zftv(ji,jj,jk+kp,jp+1) = zftv(ji,jj,jk+kp,jp+1) + & 
     219                     &                     zah * r1_e2v(ji,jj-jp) * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 
    194220                  ! 
    195221               END_3D 
     
    197223         END DO 
    198224         ! 
     225         ! NOTE: [halo1-halo2] Use DO_3D instead of DO_3D_OVR 
     226         ! jp loop contributions added here to ensure consistent floating point arithmetic for different nn_hls 
     227         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
     228            ah_wslp2(ji,jj,jk) = ah_wslp2(ji,jj,jk) + (zftu(ji,jj,jk,1) + zftu(ji,jj,jk,2)) 
     229            akz     (ji,jj,jk) = akz     (ji,jj,jk) + (zftv(ji,jj,jk,1) + zftv(ji,jj,jk,2)) 
     230         END_3D 
     231         ! 
    199232         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient 
    200233            ! 
    201234            IF( ln_traldf_blp ) THEN                ! bilaplacian operator 
    202235               ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    203                DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  
     236               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    204237                  akz(ji,jj,jk) = 16._wp           & 
    205238                     &   * ah_wslp2   (ji,jj,jk)   & 
     
    210243            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator 
    211244               ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
    212                DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  
     245               DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 
    213246                  ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 
    214247                  zcoef0 = rDt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  ) 
     
    219252         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 
    220253            ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpk ) 
    221             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )  
     254            DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 
    222255               akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 
    223256            END_3D 
    224257         ENDIF 
    225258         ! 
    226          ! TEMP: [tiling] These changes not necessary if XIOS has subdomain support 
    227          IF( .NOT. l_istiled .OR. ntile == nijtile )  THEN                ! Do only for the full domain 
    228             IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 
    229                IF( ln_tile ) CALL dom_tile_stop( ldhold=.TRUE., cstr='traldf_triad' ) 
    230  
    231                zpsi_uw(:,:,:) = 0._wp 
    232                zpsi_vw(:,:,:) = 0._wp 
    233  
    234                DO jp = 0, 1 
    235                   DO kp = 0, 1 
    236                      ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
    237                      DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )  
    238                         zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 
    239                            & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp) 
    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+jp,jk,1-jp,kp) 
    242                      END_3D 
    243                   END DO 
     259         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 
     260            zpsi_uw(:,:,:) = 0._wp 
     261            zpsi_vw(:,:,:) = 0._wp 
     262 
     263            DO jp = 0, 1 
     264               DO kp = 0, 1 
     265                  DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 
     266                     zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 
     267                        & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp) 
     268                     zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 
     269                        & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+jp,jk,1-jp,kp) 
     270                  END_3D 
    244271               END DO 
    245                CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 
    246  
    247                IF( ln_tile .AND. .NOT. l_istiled ) CALL dom_tile_start( ldhold=.TRUE., cstr='traldf_triad' ) 
    248             ENDIF 
     272            END DO 
     273            CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 
    249274         ENDIF 
    250275         ! 
     
    256281         ! Zero fluxes for each tracer 
    257282!!gm  this should probably be done outside the jn loop 
    258          ztfw(:,:,:) = 0._wp 
    259          zftu(:,:,:) = 0._wp 
    260          zftv(:,:,:) = 0._wp 
     283         ! NOTE: [halo1-halo2] 
     284         ztfw(:,:,:,:,:) = 0._wp 
     285         zftu(:,:,:,:) = 0._wp 
     286         zftv(:,:,:,:) = 0._wp 
    261287         ! 
    262288         ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 )    !==  before lateral T & S gradients at T-level jk  ==! 
    263          DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )  
     289         ! NOTE: [halo1-halo2] 
     290         DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) 
    264291            zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 
    265292            zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 
     
    267294         IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level 
    268295            ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 )                    ! bottom level 
    269             DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )  
     296            ! NOTE: [halo1-halo2] 
     297            DO_2D( iij, iij-1, iij, iij-1 ) 
    270298               zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 
    271299               zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     
    273301            IF( ln_isfcav ) THEN                   ! top level (ocean cavities only) 
    274302               ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
    275                DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )  
     303               ! NOTE: [halo1-halo2] 
     304               DO_2D( iij, iij-1, iij, iij-1 ) 
    276305                  IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 
    277306                  IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 
     
    287316            !                    !==  Vertical tracer gradient at level jk and jk+1 
    288317            ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) 
    289             DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )  
     318            ! NOTE: [halo1-halo2] 
     319            DO_2D( iij, iij, iij, iij ) 
    290320               zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 
    291321            END_2D 
     
    295325            ELSE 
    296326               ! [comm_cleanup] ! DO_2D( 1, 1, 1, 1 ) 
    297                DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )  
     327               ! NOTE: [halo1-halo2] 
     328               DO_2D( iij, iij, iij, iij ) 
    298329                  zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    299330               END_2D 
     
    306337                  DO kp = 0, 1 
    307338                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
    308                      DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )  
     339                     ! NOTE: [halo1-halo2] 
     340                     DO_2D( iij, iij-1, iij, iij-1 ) 
    309341                        ze1ur = r1_e1u(ji,jj) 
    310342                        zdxt  = zdit(ji,jj,jk) * ze1ur 
     
    319351                        zah_slp  = zah * zslope_iso 
    320352                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew 
    321                         zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
    322                         ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt                 * zbu * ze3wr 
     353                        ! NOTE: [halo1-halo2] 
     354                        zftu(ji   ,jj,jk,ip+1  ) = zftu(ji   ,jj,jk,ip+1  ) - & 
     355                           &                       ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
     356                        ztfw(ji,jj,jk+kp,ip+1,1) = ztfw(ji,jj,jk+kp,ip+1,1) - & 
     357                           &                                      (zah_slp + zaei_slp) * zdxt   * zbu * ze3wr 
    323358                     END_2D 
    324359                  END DO 
     
    328363                  DO kp = 0, 1 
    329364                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
    330                      DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )  
     365                     ! NOTE: [halo1-halo2] 
     366                     DO_2D( iij, iij-1, iij, iij-1 ) 
    331367                        ze2vr = r1_e2v(ji,jj) 
    332368                        zdyt  = zdjt(ji,jj,jk) * ze2vr 
     
    340376                        zah_slp = zah * zslope_iso 
    341377                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew 
    342                         zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
    343                         ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt                * zbv * ze3wr 
     378                        ! NOTE: [halo1-halo2] 
     379                        zftv(ji,jj   ,jk,jp+1  ) = zftv(ji,jj   ,jk,jp+1  ) - & 
     380                           &                       ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
     381                        ztfw(ji,jj,jk+kp,jp+1,2) = ztfw(ji,jj,jk+kp,jp+1,2) - & 
     382                           &                                      (zah_slp + zaei_slp) * zdyt   * zbv * ze3wr 
    344383                     END_2D 
    345384                  END DO 
     
    351390                  DO kp = 0, 1 
    352391                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
    353                      DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )  
     392                     ! NOTE: [halo1-halo2] 
     393                     DO_2D( iij, iij-1, iij, iij-1 ) 
    354394                        ze1ur = r1_e1u(ji,jj) 
    355395                        zdxt  = zdit(ji,jj,jk) * ze1ur 
     
    364404                        zah_slp  = zah * zslope_iso 
    365405                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! aeit(ji+ip,jj,jk)*zslope_skew 
    366                         zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
    367                         ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 
     406                        ! NOTE: [halo1-halo2] 
     407                        zftu(ji   ,jj,jk,ip+1  ) = zftu(ji   ,jj,jk,ip+1  ) - & 
     408                           &                       ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
     409                        ztfw(ji,jj,jk+kp,ip+1,1) = ztfw(ji,jj,jk+kp,ip+1,1) - & 
     410                           &                                      (zah_slp + zaei_slp) * zdxt   * zbu * ze3wr 
    368411                     END_2D 
    369412                  END DO 
     
    373416                  DO kp = 0, 1 
    374417                     ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 
    375                      DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )  
     418                     ! NOTE: [halo1-halo2] 
     419                     DO_2D( iij, iij-1, iij, iij-1 ) 
    376420                        ze2vr = r1_e2v(ji,jj) 
    377421                        zdyt  = zdjt(ji,jj,jk) * ze2vr 
     
    385429                        zah_slp = zah * zslope_iso 
    386430                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! aeit(ji,jj+jp,jk)*zslope_skew 
    387                         zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
    388                         ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 
     431                        ! NOTE: [halo1-halo2] 
     432                        zftv(ji,jj,jk,   jp+1  ) = zftv(ji,jj,jk,   jp+1  ) - & 
     433                           &                       ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
     434                        ztfw(ji,jj,jk+kp,jp+1,2) = ztfw(ji,jj,jk+kp,jp+1,2) - & 
     435                           &                                      (zah_slp + zaei_slp) * zdyt   * zbv * ze3wr 
    389436                     END_2D 
    390437                  END DO 
    391438               END DO 
    392439            ENDIF 
     440            ! NOTE: [halo1-halo2] 
     441            ! ip and jp loop contributions added here to ensure consistent floating point arithmetic for different nn_hls 
     442            DO_2D( iij, iij-1, iij, iij-1 ) 
     443               zftu(ji,jj,jk,1) = zftu(ji,jj,jk,1) + zftu(ji,jj,jk,2) 
     444               zftv(ji,jj,jk,1) = zftv(ji,jj,jk,1) + zftv(ji,jj,jk,2) 
     445            END_2D 
     446            DO_2D( iij-1, iij-1, iij-1, iij-1 ) 
     447               ztfw(ji,jj,jk,1,1) = (ztfw(ji,jj,jk,1,1) + ztfw(ji-1,jj,jk,2,1)) + & 
     448                  &                 (ztfw(ji,jj,jk,1,2) + ztfw(ji,jj-1,jk,2,2)) 
     449            END_2D 
    393450            !                             !==  horizontal divergence and add to the general trend  ==! 
    394451            ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 
    395             DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )  
     452            ! NOTE: [halo1-halo2] 
     453            DO_2D( iij-1, iij-1, iij-1, iij-1 ) 
     454               ! NOTE: [halo1-halo2] 
     455               ! Extra brackets required to ensure consistent floating point arithmetic for different nn_hls for bilaplacian 
    396456               pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    & 
    397                   &                       + zsign * (  zftu(ji-1,jj  ,jk) - zftu(ji,jj,jk)       & 
    398                   &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   & 
    399                   &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  ) 
     457                  &                       + zsign * (  (zftu(ji-1,jj,jk,1) - zftu(ji,jj,jk,1))       & 
     458                  &                                  + (zftv(ji,jj-1,jk,1) - zftv(ji,jj,jk,1))   )   & 
     459                  &                               / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  ) 
    400460            END_2D 
    401461            ! 
     
    405465         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
    406466            ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 
    407             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  
    408                ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   & 
     467            ! NOTE: [halo1-halo2] 
     468            DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 
     469               ztfw(ji,jj,jk,1,1) = ztfw(ji,jj,jk,1,1) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   & 
    409470                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    410471                  &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
     
    414475            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    415476               ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 
    416                DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  
    417                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             & 
     477               ! NOTE: [halo1-halo2] 
     478               DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 
     479                  ztfw(ji,jj,jk,1,1) = ztfw(ji,jj,jk,1,1) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             & 
    418480                     &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    419481               END_3D 
    420482            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp. 
    421                ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 
    422                DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )  
    423                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      & 
     483               ! NOTE: [halo1-halo2] 
     484               DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
     485                  ztfw(ji,jj,jk,1,1) = ztfw(ji,jj,jk,1,1) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      & 
    424486                     &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   & 
    425487                     &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   ) 
     
    429491         ! 
    430492         ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 )      !==  Divergence of vertical fluxes added to pta  ==! 
    431          DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )  
    432             pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn)    & 
    433             &                                  + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   & 
     493         ! NOTE: [halo1-halo2] 
     494         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) 
     495            pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1,1,1) - ztfw(ji,jj,jk,1,1)  )   & 
    434496               &                                              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 
    435497         END_3D 
     
    439501            ! 
    440502            !                          ! "Poleward" diffusive heat or salt transports (T-S case only) 
    441             IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:)  ) 
     503            ! NOTE: [halo1-halo2] 
     504            IF( l_ptr )  CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:,1)  ) 
    442505            !                          ! Diffusive heat transports 
    443             IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) ) 
     506            ! NOTE: [halo1-halo2] 
     507            IF( l_hst )  CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:,1), zftv(:,:,:,1) ) 
    444508            ! 
    445509         ENDIF                                                    !== end pass selection  ==! 
Note: See TracChangeset for help on using the changeset viewer.