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 12340 for NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90 – NEMO

Ignore:
Timestamp:
2020-01-27T15:31:53+01:00 (4 years ago)
Author:
acc
Message:

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/icedyn_adv_umx.F90

    r12252 r12340  
    5252   !! * Substitutions 
    5353#  include "vectopt_loop_substitute.h90" 
     54#  include "do_loop_substitute.h90" 
    5455   !!---------------------------------------------------------------------- 
    5556   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    … …  
    107108      ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
    108109      DO jl = 1, jpl 
    109          DO jj = 2, jpjm1 
    110             DO ji = fs_2, fs_jpim1 
    111                zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
    112                   &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
    113                   &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
    114                   &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
    115                zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
    116                   &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
    117                   &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
    118                   &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
    119                zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
    120                   &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
    121                   &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    122                   &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    123             END DO 
    124          END DO 
     110         DO_2D_00_00 
     111            zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), & 
     112               &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), & 
     113               &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 
     114               &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 
     115            zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), & 
     116               &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), & 
     117               &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 
     118               &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 
     119            zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), & 
     120               &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), & 
     121               &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
     122               &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
     123         END_2D 
    125124      END DO 
    126125      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1. ) 
    … …  
    152151      ! 
    153152      ! --- define velocity for advection: u*grad(H) --- ! 
    154       DO jj = 2, jpjm1 
    155          DO ji = fs_2, fs_jpim1 
    156             IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp 
    157             ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj) 
    158             ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj) 
    159             ENDIF 
    160  
    161             IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp 
    162             ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1) 
    163             ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  ) 
    164             ENDIF 
    165          END DO 
    166       END DO 
     153      DO_2D_00_00 
     154         IF    ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN   ;   zcu_box(ji,jj) = 0._wp 
     155         ELSEIF( pu_ice(ji,jj)                   >  0._wp ) THEN   ;   zcu_box(ji,jj) = pu_ice(ji-1,jj) 
     156         ELSE                                                      ;   zcu_box(ji,jj) = pu_ice(ji  ,jj) 
     157         ENDIF 
     158 
     159         IF    ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN   ;   zcv_box(ji,jj) = 0._wp 
     160         ELSEIF( pv_ice(ji,jj)                   >  0._wp ) THEN   ;   zcv_box(ji,jj) = pv_ice(ji,jj-1) 
     161         ELSE                                                      ;   zcv_box(ji,jj) = pv_ice(ji,jj  ) 
     162         ENDIF 
     163      END_2D 
    167164 
    168165      !---------------! 
    … …  
    187184            IF( .NOT. ALLOCATED(jmsk_small) )   ALLOCATE( jmsk_small(jpi,jpj,jpl) )  
    188185            DO jl = 1, jpl 
    189                DO jj = 1, jpjm1 
    190                   DO ji = 1, jpim1 
    191                      zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) ) 
    192                      IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0 
    193                      ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF 
    194                      zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) ) 
    195                      IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0 
    196                      ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF 
    197                   END DO 
    198                END DO 
     186               DO_2D_10_10 
     187                  zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) ) 
     188                  IF( zvi_cen < epsi06) THEN   ;   imsk_small(ji,jj,jl) = 0 
     189                  ELSE                         ;   imsk_small(ji,jj,jl) = 1   ;   ENDIF 
     190                  zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) ) 
     191                  IF( zvi_cen < epsi06) THEN   ;   jmsk_small(ji,jj,jl) = 0 
     192                  ELSE                         ;   jmsk_small(ji,jj,jl) = 1   ;   ENDIF 
     193               END_2D 
    199194            END DO 
    200195         ENDIF 
    … …  
    338333         !== Open water area ==! 
    339334         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 
    340          DO jj = 2, jpjm1 
    341             DO ji = fs_2, fs_jpim1 
    342                pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &  
    343                   &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
    344             END DO 
    345          END DO 
     335         DO_2D_00_00 
     336            pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &  
     337               &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 
     338         END_2D 
    346339         CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T',  1. ) 
    347340         ! 
    … …  
    449442      IF( pamsk == 0._wp ) THEN 
    450443         DO jl = 1, jpl 
    451             DO jj = 1, jpjm1 
    452                DO ji = 1, fs_jpim1 
    453                   IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 
    454                      zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj) 
    455                      zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 
    456                   ELSE 
    457                      zfu_ho (ji,jj,jl) = 0._wp 
    458                      zfu_ups(ji,jj,jl) = 0._wp 
    459                   ENDIF 
    460                   ! 
    461                   IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 
    462                      zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj) 
    463                      zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 
    464                   ELSE 
    465                      zfv_ho (ji,jj,jl) = 0._wp   
    466                      zfv_ups(ji,jj,jl) = 0._wp   
    467                   ENDIF 
    468                END DO 
    469             END DO 
     444            DO_2D_10_10 
     445               IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 
     446                  zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj) 
     447                  zfu_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) * pua_ups(ji,jj,jl) / pu(ji,jj) 
     448               ELSE 
     449                  zfu_ho (ji,jj,jl) = 0._wp 
     450                  zfu_ups(ji,jj,jl) = 0._wp 
     451               ENDIF 
     452               ! 
     453               IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 
     454                  zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj) 
     455                  zfv_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) * pva_ups(ji,jj,jl) / pv(ji,jj) 
     456               ELSE 
     457                  zfv_ho (ji,jj,jl) = 0._wp   
     458                  zfv_ups(ji,jj,jl) = 0._wp   
     459               ENDIF 
     460            END_2D 
    470461         END DO 
    471462 
    … …  
    473464         ! thus we calculate the upstream solution and apply a limiter again 
    474465         DO jl = 1, jpl 
    475             DO jj = 2, jpjm1 
    476                DO ji = fs_2, fs_jpim1 
    477                   ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 
    478                   ! 
    479                   zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 
    480                END DO 
    481             END DO 
     466            DO_2D_00_00 
     467               ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) + zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) 
     468               ! 
     469               zt_ups(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 
     470            END_2D 
    482471         END DO 
    483472         CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T',  1. ) 
    … …  
    496485      IF( PRESENT( pua_ho ) ) THEN 
    497486         DO jl = 1, jpl 
    498             DO jj = 1, jpjm1 
    499                DO ji = 1, fs_jpim1 
    500                   pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
    501                   pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
    502               END DO 
    503             END DO 
     487            DO_2D_10_10 
     488               pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 
     489               pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 
     490            END_2D 
    504491         END DO 
    505492      ENDIF 
    … …  
    508495      ! --------------------------------- 
    509496      DO jl = 1, jpl 
    510          DO jj = 2, jpjm1 
    511             DO ji = fs_2, fs_jpim1  
    512                ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )   
    513                ! 
    514                ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)                
    515             END DO 
    516          END DO 
     497         DO_2D_00_00 
     498            ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) )   
     499            ! 
     500            ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1)                
     501         END_2D 
    517502      END DO 
    518503      CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1. ) 
    … …  
    544529         ! 
    545530         DO jl = 1, jpl 
    546             DO jj = 1, jpjm1 
    547                DO ji = 1, fs_jpim1 
     531            DO_2D_10_10 
     532               pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     533               pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
     534            END_2D 
     535         END DO 
     536         ! 
     537      ELSE                              !** alternate directions **! 
     538         ! 
     539         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     540            ! 
     541            DO jl = 1, jpl              !-- flux in x-direction 
     542               DO_2D_10_10 
    548543                  pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
     544               END_2D 
     545            END DO 
     546            ! 
     547            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
     548               DO_2D_00_00 
     549                  ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              & 
     550                     &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     551                  ! 
     552                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     553               END_2D 
     554            END DO 
     555            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
     556            ! 
     557            DO jl = 1, jpl              !-- flux in y-direction 
     558               DO_2D_10_10 
     559                  pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 
     560               END_2D 
     561            END DO 
     562            ! 
     563         ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
     564            ! 
     565            DO jl = 1, jpl              !-- flux in y-direction 
     566               DO_2D_10_10 
    549567                  pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
    550                END DO 
    551             END DO 
    552          END DO 
    553          ! 
    554       ELSE                              !** alternate directions **! 
    555          ! 
    556          IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==! 
     568               END_2D 
     569            END DO 
     570            ! 
     571            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
     572               DO_2D_00_00 
     573                  ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  & 
     574                     &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     575                  ! 
     576                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     577               END_2D 
     578            END DO 
     579            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    557580            ! 
    558581            DO jl = 1, jpl              !-- flux in x-direction 
    559                DO jj = 1, jpjm1 
    560                   DO ji = 1, fs_jpim1 
    561                      pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 
    562                   END DO 
    563                END DO 
    564             END DO 
    565             ! 
    566             DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    567                DO jj = 2, jpjm1 
    568                   DO ji = fs_2, fs_jpim1 
    569                      ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              & 
    570                         &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    571                      ! 
    572                      zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
    573                   END DO 
    574                END DO 
    575             END DO 
    576             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    577             ! 
    578             DO jl = 1, jpl              !-- flux in y-direction 
    579                DO jj = 1, jpjm1 
    580                   DO ji = 1, fs_jpim1 
    581                      pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 
    582                   END DO 
    583                END DO 
    584             END DO 
    585             ! 
    586          ELSE                                                               !==  even ice time step:  adv_y then adv_x  ==! 
    587             ! 
    588             DO jl = 1, jpl              !-- flux in y-direction 
    589                DO jj = 1, jpjm1 
    590                   DO ji = 1, fs_jpim1 
    591                      pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 
    592                   END DO 
    593                END DO 
    594             END DO 
    595             ! 
    596             DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    597                DO jj = 2, jpjm1 
    598                   DO ji = fs_2, fs_jpim1 
    599                      ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  & 
    600                         &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    601                      ! 
    602                      zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
    603                   END DO 
    604                END DO 
    605             END DO 
    606             CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    607             ! 
    608             DO jl = 1, jpl              !-- flux in x-direction 
    609                DO jj = 1, jpjm1 
    610                   DO ji = 1, fs_jpim1 
    611                      pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 
    612                   END DO 
    613                END DO 
     582               DO_2D_10_10 
     583                  pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 
     584               END_2D 
    614585            END DO 
    615586            ! 
    … …  
    619590      ! 
    620591      DO jl = 1, jpl                    !-- after tracer with upstream scheme 
    621          DO jj = 2, jpjm1 
    622             DO ji = fs_2, fs_jpim1 
    623                ztra = - (   pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj  ,jl)   & 
    624                   &       + pfv_ups(ji,jj,jl) - pfv_ups(ji  ,jj-1,jl) ) & 
    625                   &   + (   pu     (ji,jj   ) - pu     (ji-1,jj     )   & 
    626                   &       + pv     (ji,jj   ) - pv     (ji  ,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    627                ! 
    628                pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
    629             END DO 
    630          END DO 
     592         DO_2D_00_00 
     593            ztra = - (   pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj  ,jl)   & 
     594               &       + pfv_ups(ji,jj,jl) - pfv_ups(ji  ,jj-1,jl) ) & 
     595               &   + (   pu     (ji,jj   ) - pu     (ji-1,jj     )   & 
     596               &       + pv     (ji,jj   ) - pv     (ji  ,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     597            ! 
     598            pt_ups(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     599         END_2D 
    631600      END DO 
    632601      CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. ) 
    … …  
    660629         ! 
    661630         DO jl = 1, jpl 
    662             DO jj = 1, jpjm1 
    663                DO ji = 1, fs_jpim1 
    664                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) ) 
    665                   pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) ) 
    666                END DO 
    667             END DO 
     631            DO_2D_10_10 
     632               pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) ) 
     633               pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) ) 
     634            END_2D 
    668635         END DO 
    669636         ! 
    … …  
    680647            ! 
    681648            DO jl = 1, jpl              !-- flux in x-direction 
    682                DO jj = 1, jpjm1 
    683                   DO ji = 1, fs_jpim1 
    684                      pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
    685                   END DO 
    686                END DO 
     649               DO_2D_10_10 
     650                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 
     651               END_2D 
    687652            END DO 
    688653            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    689654 
    690655            DO jl = 1, jpl              !-- first guess of tracer from u-flux 
    691                DO jj = 2, jpjm1 
    692                   DO ji = fs_2, fs_jpim1 
    693                      ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              & 
    694                         &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    695                      ! 
    696                      zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
    697                   END DO 
    698                END DO 
     656               DO_2D_00_00 
     657                  ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              & 
     658                     &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     659                  ! 
     660                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     661               END_2D 
    699662            END DO 
    700663            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    701664 
    702665            DO jl = 1, jpl              !-- flux in y-direction 
    703                DO jj = 1, jpjm1 
    704                   DO ji = 1, fs_jpim1 
    705                      pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 
    706                   END DO 
    707                END DO 
     666               DO_2D_10_10 
     667                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 
     668               END_2D 
    708669            END DO 
    709670            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    … …  
    712673            ! 
    713674            DO jl = 1, jpl              !-- flux in y-direction 
    714                DO jj = 1, jpjm1 
    715                   DO ji = 1, fs_jpim1 
    716                      pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
    717                   END DO 
    718                END DO 
     675               DO_2D_10_10 
     676                  pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 
     677               END_2D 
    719678            END DO 
    720679            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 
    721680            ! 
    722681            DO jl = 1, jpl              !-- first guess of tracer from v-flux 
    723                DO jj = 2, jpjm1 
    724                   DO ji = fs_2, fs_jpim1 
    725                      ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  & 
    726                         &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
    727                      ! 
    728                      zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
    729                   END DO 
    730                END DO 
     682               DO_2D_00_00 
     683                  ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  & 
     684                     &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk) 
     685                  ! 
     686                  zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     687               END_2D 
    731688            END DO 
    732689            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    733690            ! 
    734691            DO jl = 1, jpl              !-- flux in x-direction 
    735                DO jj = 1, jpjm1 
    736                   DO ji = 1, fs_jpim1 
    737                      pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 
    738                   END DO 
    739                END DO 
     692               DO_2D_10_10 
     693                  pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 
     694               END_2D 
    740695            END DO 
    741696            IF( np_limiter == 2 .OR. np_limiter == 3 )   CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 
    … …  
    783738         !                                                        !--  advective form update in zpt  --! 
    784739         DO jl = 1, jpl 
    785             DO jj = 2, jpjm1 
    786                DO ji = fs_2, fs_jpim1 
    787                   zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pubox(ji,jj   ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t  (ji,jj) & 
    788                      &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) & 
    789                      &                                                                                        * pamsk           & 
    790                      &                             ) * pdt ) * tmask(ji,jj,1) 
    791                END DO 
    792             END DO 
     740            DO_2D_00_00 
     741               zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pubox(ji,jj   ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t  (ji,jj) & 
     742                  &                              + pt   (ji,jj,jl) * ( pu  (ji,jj   ) - pu  (ji-1,jj   ) ) * r1_e1e2t(ji,jj) & 
     743                  &                                                                                        * pamsk           & 
     744                  &                             ) * pdt ) * tmask(ji,jj,1) 
     745            END_2D 
    793746         END DO 
    794747         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    … …  
    812765         !                                                        !--  advective form update in zpt  --! 
    813766         DO jl = 1, jpl 
    814             DO jj = 2, jpjm1 
    815                DO ji = fs_2, fs_jpim1 
    816                   zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pvbox(ji,jj   ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t  (ji,jj) & 
    817                      &                              + pt   (ji,jj,jl) * ( pv  (ji,jj   ) - pv  (ji,jj-1   ) ) * r1_e1e2t(ji,jj) & 
    818                      &                                                                                        * pamsk           & 
    819                      &                             ) * pdt ) * tmask(ji,jj,1)  
    820                END DO 
    821             END DO 
     767            DO_2D_00_00 
     768               zpt(ji,jj,jl) = ( pt(ji,jj,jl) - (  pvbox(ji,jj   ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t  (ji,jj) & 
     769                  &                              + pt   (ji,jj,jl) * ( pv  (ji,jj   ) - pv  (ji,jj-1   ) ) * r1_e1e2t(ji,jj) & 
     770                  &                                                                                        * pamsk           & 
     771                  &                             ) * pdt ) * tmask(ji,jj,1)  
     772            END_2D 
    822773         END DO 
    823774         CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 
    … …  
    896847         !         
    897848         DO jl = 1, jpl 
    898             DO jj = 1, jpjm1 
    899                DO ji = 1, fs_jpim1   ! vector opt. 
    900                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    901                      &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
    902                END DO 
    903             END DO 
     849            DO_2D_10_10 
     850               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     851                  &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     852            END_2D 
    904853         END DO 
    905854         ! 
    … …  
    907856         ! 
    908857         DO jl = 1, jpl 
    909             DO jj = 1, jpjm1 
    910                DO ji = 1, fs_jpim1   ! vector opt. 
    911                   zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    912                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    913                      &                                                            - zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )  
    914                END DO 
    915             END DO 
     858            DO_2D_10_10 
     859               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     860               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     861                  &                                                            - zcu   * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )  
     862            END_2D 
    916863         END DO 
    917864         !   
    … …  
    919866         ! 
    920867         DO jl = 1, jpl 
    921             DO jj = 1, jpjm1 
    922                DO ji = 1, fs_jpim1   ! vector opt. 
    923                   zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    924                   zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     868            DO_2D_10_10 
     869               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     870               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    925871!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    926                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     & 
    927                      &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) & 
    928                      &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     & 
    929                      &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 
    930                END DO 
    931             END DO 
     872               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     & 
     873                  &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) & 
     874                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     & 
     875                  &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 
     876            END_2D 
    932877         END DO 
    933878         ! 
    … …  
    935880         ! 
    936881         DO jl = 1, jpl 
    937             DO jj = 1, jpjm1 
    938                DO ji = 1, fs_jpim1   ! vector opt. 
    939                   zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    940                   zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     882            DO_2D_10_10 
     883               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     884               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    941885!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    942                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     & 
    943                      &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) & 
    944                      &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     & 
    945                      &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 
    946                END DO 
    947             END DO 
     886               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     & 
     887                  &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) & 
     888                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     & 
     889                  &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 
     890            END_2D 
    948891         END DO 
    949892         ! 
    … …  
    951894         ! 
    952895         DO jl = 1, jpl 
    953             DO jj = 1, jpjm1 
    954                DO ji = 1, fs_jpim1   ! vector opt. 
    955                   zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
    956                   zdx2 = e1u(ji,jj) * e1u(ji,jj) 
     896            DO_2D_10_10 
     897               zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 
     898               zdx2 = e1u(ji,jj) * e1u(ji,jj) 
    957899!!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj) 
    958                   zdx4 = zdx2 * zdx2 
    959                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (        (                       pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     & 
    960                      &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) & 
    961                      &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) * (                       ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     & 
    962                      &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 
    963                      &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)     & 
    964                      &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 
    965                END DO 
    966             END DO 
     900               zdx4 = zdx2 * zdx2 
     901               pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (        (                       pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     & 
     902                  &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) & 
     903                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) * (                       ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     & 
     904                  &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 
     905                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)     & 
     906                  &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 
     907            END_2D 
    967908         END DO 
    968909         ! 
    … …  
    974915      IF( ll_neg ) THEN 
    975916         DO jl = 1, jpl 
    976             DO jj = 1, jpjm1 
    977                DO ji = 1, fs_jpim1 
    978                   IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    979                      pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
    980                         &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
    981                   ENDIF 
    982                END DO 
    983             END DO 
     917            DO_2D_10_10 
     918               IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
     919                  pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   & 
     920                     &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 
     921               ENDIF 
     922            END_2D 
    984923         END DO 
    985924      ENDIF 
    986925      !                                                     !-- High order flux in i-direction  --! 
    987926      DO jl = 1, jpl 
    988          DO jj = 1, jpjm1 
    989             DO ji = 1, fs_jpim1   ! vector opt. 
    990                pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 
    991             END DO 
    992          END DO 
     927         DO_2D_10_10 
     928            pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl) 
     929         END_2D 
    993930      END DO 
    994931      ! 
    … …  
    1021958      !                                                     !--  Laplacian in j-direction  --! 
    1022959      DO jl = 1, jpl 
    1023          DO jj = 1, jpjm1         ! First derivative (gradient) 
    1024             DO ji = fs_2, fs_jpim1 
    1025                ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    1026             END DO 
    1027          END DO 
    1028          DO jj = 2, jpjm1         ! Second derivative (Laplacian) 
    1029             DO ji = fs_2, fs_jpim1 
    1030                ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
    1031             END DO 
    1032          END DO 
     960         DO_2D_10_00 
     961            ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     962         END_2D 
     963         DO_2D_00_00 
     964            ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
     965         END_2D 
    1033966      END DO 
    1034967      CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1. ) 
    … …  
    1036969      !                                                     !--  BiLaplacian in j-direction  --! 
    1037970      DO jl = 1, jpl 
    1038          DO jj = 1, jpjm1         ! First derivative 
    1039             DO ji = fs_2, fs_jpim1 
    1040                ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    1041             END DO 
    1042          END DO 
    1043          DO jj = 2, jpjm1         ! Second derivative 
    1044             DO ji = fs_2, fs_jpim1 
    1045                ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
    1046             END DO 
    1047          END DO 
     971         DO_2D_10_00 
     972            ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
     973         END_2D 
     974         DO_2D_00_00 
     975            ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
     976         END_2D 
    1048977      END DO 
    1049978      CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1. ) 
    … …  
    1054983      CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21) 
    1055984         DO jl = 1, jpl 
    1056             DO jj = 1, jpjm1 
    1057                DO ji = 1, fs_jpim1 
    1058                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
    1059                      &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    1060                END DO 
    1061             END DO 
     985            DO_2D_10_10 
     986               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
     987                  &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     988            END_2D 
    1062989         END DO 
    1063990         ! 
    1064991      CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23) 
    1065992         DO jl = 1, jpl 
    1066             DO jj = 1, jpjm1 
    1067                DO ji = 1, fs_jpim1 
    1068                   zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    1069                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
    1070                      &                                                            - zcv *   ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    1071                END DO 
    1072             END DO 
     993            DO_2D_10_10 
     994               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     995               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   & 
     996                  &                                                            - zcv *   ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     997            END_2D 
    1073998         END DO 
    1074999         ! 
    10751000      CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24) 
    10761001         DO jl = 1, jpl 
    1077             DO jj = 1, jpjm1 
    1078                DO ji = 1, fs_jpim1 
    1079                   zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    1080                   zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     1002            DO_2D_10_10 
     1003               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1004               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    10811005!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    1082                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     & 
    1083                      &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) & 
    1084                      &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     & 
    1085                      &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 
    1086                END DO 
    1087             END DO 
     1006               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     & 
     1007                  &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) & 
     1008                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     & 
     1009                  &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 
     1010            END_2D 
    10881011         END DO 
    10891012         ! 
    10901013      CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27) 
    10911014         DO jl = 1, jpl 
    1092             DO jj = 1, jpjm1 
    1093                DO ji = 1, fs_jpim1 
    1094                   zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    1095                   zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     1015            DO_2D_10_10 
     1016               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1017               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    10961018!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    1097                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     & 
    1098                      &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) & 
    1099                      &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     & 
    1100                      &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 
    1101                END DO 
    1102             END DO 
     1019               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     & 
     1020                  &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) & 
     1021                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     & 
     1022                  &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 
     1023            END_2D 
    11031024         END DO 
    11041025         ! 
    11051026      CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29) 
    11061027         DO jl = 1, jpl 
    1107             DO jj = 1, jpjm1 
    1108                DO ji = 1, fs_jpim1 
    1109                   zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
    1110                   zdy2 = e2v(ji,jj) * e2v(ji,jj) 
     1028            DO_2D_10_10 
     1029               zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 
     1030               zdy2 = e2v(ji,jj) * e2v(ji,jj) 
    11111031!!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj) 
    1112                   zdy4 = zdy2 * zdy2 
    1113                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     & 
    1114                      &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) & 
    1115                      &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) * (                       ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     & 
    1116                      &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 
    1117                      &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     & 
    1118                      &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 
    1119                END DO 
    1120             END DO 
     1032               zdy4 = zdy2 * zdy2 
     1033               pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     & 
     1034                  &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) & 
     1035                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) * (                       ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     & 
     1036                  &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 
     1037                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     & 
     1038                  &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 
     1039            END_2D 
    11211040         END DO 
    11221041         ! 
    … …  
    11281047      IF( ll_neg ) THEN 
    11291048         DO jl = 1, jpl 
    1130             DO jj = 1, jpjm1 
    1131                DO ji = 1, fs_jpim1 
    1132                   IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
    1133                      pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
    1134                         &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
    1135                   ENDIF 
    1136                END DO 
    1137             END DO 
     1049            DO_2D_10_10 
     1050               IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 
     1051                  pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  & 
     1052                     &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 
     1053               ENDIF 
     1054            END_2D 
    11381055         END DO 
    11391056      ENDIF 
    11401057      !                                                     !-- High order flux in j-direction  --! 
    11411058      DO jl = 1, jpl 
    1142          DO jj = 1, jpjm1 
    1143             DO ji = 1, fs_jpim1   ! vector opt. 
    1144                pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 
    1145             END DO 
    1146          END DO 
     1059         DO_2D_10_10 
     1060            pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 
     1061         END_2D 
    11471062      END DO 
    11481063      ! 
    … …  
    11781093      ! -------------------------------------------------- 
    11791094      DO jl = 1, jpl 
    1180          DO jj = 1, jpjm1 
    1181             DO ji = 1, fs_jpim1   ! vector opt. 
    1182                pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 
    1183                pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 
    1184             END DO 
    1185          END DO 
     1095         DO_2D_10_10 
     1096            pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 
     1097            pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 
     1098         END_2D 
    11861099      END DO 
    11871100 
    … …  
    11971110          
    11981111         DO jl = 1, jpl 
    1199             DO jj = 2, jpjm1 
    1200                DO ji = fs_2, fs_jpim1  
    1201                   zti_ups(ji,jj,jl)= pt_ups(ji+1,jj  ,jl) 
    1202                   ztj_ups(ji,jj,jl)= pt_ups(ji  ,jj+1,jl) 
    1203                END DO 
    1204             END DO 
     1112            DO_2D_00_00 
     1113               zti_ups(ji,jj,jl)= pt_ups(ji+1,jj  ,jl) 
     1114               ztj_ups(ji,jj,jl)= pt_ups(ji  ,jj+1,jl) 
     1115            END_2D 
    12051116         END DO 
    12061117         CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. ) 
    12071118 
    12081119         DO jl = 1, jpl 
    1209             DO jj = 2, jpjm1 
    1210                DO ji = fs_2, fs_jpim1 
    1211                   IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj  ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND.  & 
    1212                      & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN 
    1213                      ! 
    1214                      IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj  ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND.  & 
    1215                         & pfv_ho(ji,jj,jl) * ( ztj_ups(ji  ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN 
    1216                         pfu_ho(ji,jj,jl)=0._wp 
    1217                         pfv_ho(ji,jj,jl)=0._wp 
    1218                      ENDIF 
    1219                      ! 
    1220                      IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj  ,jl) ) <= 0._wp .AND.  & 
    1221                         & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji  ,jj-1,jl) ) <= 0._wp ) THEN 
    1222                         pfu_ho(ji,jj,jl)=0._wp 
    1223                         pfv_ho(ji,jj,jl)=0._wp 
    1224                      ENDIF 
    1225                      ! 
     1120            DO_2D_00_00 
     1121               IF ( pfu_ho(ji,jj,jl) * ( pt_ups(ji+1,jj  ,jl) - pt_ups(ji,jj,jl) ) <= 0._wp .AND.  & 
     1122                  & pfv_ho(ji,jj,jl) * ( pt_ups(ji  ,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0._wp ) THEN 
     1123                  ! 
     1124                  IF(  pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj  ,jl) - zti_ups(ji,jj,jl) ) <= 0._wp .AND.  & 
     1125                     & pfv_ho(ji,jj,jl) * ( ztj_ups(ji  ,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0._wp ) THEN 
     1126                     pfu_ho(ji,jj,jl)=0._wp 
     1127                     pfv_ho(ji,jj,jl)=0._wp 
    12261128                  ENDIF 
    1227                END DO 
    1228             END DO 
     1129                  ! 
     1130                  IF(  pfu_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji-1,jj  ,jl) ) <= 0._wp .AND.  & 
     1131                     & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj,jl) - pt_ups(ji  ,jj-1,jl) ) <= 0._wp ) THEN 
     1132                     pfu_ho(ji,jj,jl)=0._wp 
     1133                     pfv_ho(ji,jj,jl)=0._wp 
     1134                  ENDIF 
     1135                  ! 
     1136               ENDIF 
     1137            END_2D 
    12291138         END DO 
    12301139         CALL lbc_lnk_multi( 'icedyn_adv_umx', pfu_ho, 'U', -1., pfv_ho, 'V', -1. )   ! lateral boundary cond. 
    … …  
    12381147      DO jl = 1, jpl 
    12391148          
    1240          DO jj = 1, jpj 
    1241             DO ji = 1, jpi 
    1242                IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 
    1243                   zbup(ji,jj) = -zbig 
    1244                   zbdo(ji,jj) =  zbig 
    1245                ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN 
    1246                   zbup(ji,jj) = pt_ups(ji,jj,jl) 
    1247                   zbdo(ji,jj) = pt_ups(ji,jj,jl) 
    1248                ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 
    1249                   zbup(ji,jj) = pt(ji,jj,jl) 
    1250                   zbdo(ji,jj) = pt(ji,jj,jl) 
    1251                ELSE 
    1252                   zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 
    1253                   zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 
    1254                ENDIF 
    1255             END DO 
    1256          END DO 
    1257  
    1258          DO jj = 2, jpjm1 
    1259             DO ji = fs_2, fs_jpim1   ! vector opt. 
    1260                ! 
    1261                zup  = MAX( zbup(ji,jj), zbup(ji-1,jj), zbup(ji+1,jj), zbup(ji,jj-1), zbup(ji,jj+1) )  ! search max/min in neighbourhood 
    1262                zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) ) 
    1263                ! 
    1264                zpos = MAX( 0._wp, pfu_ho(ji-1,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji  ,jj  ,jl) ) &  ! positive/negative part of the flux 
    1265                   & + MAX( 0._wp, pfv_ho(ji  ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj  ,jl) ) 
    1266                zneg = MAX( 0._wp, pfu_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj  ,jl) ) & 
    1267                   & + MAX( 0._wp, pfv_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj-1,jl) ) 
    1268                ! 
    1269                zpos = zpos - (pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
    1270                   &          ) * ( 1. - pamsk ) 
    1271                zneg = zneg + (pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
    1272                   &          ) * ( 1. - pamsk ) 
    1273                ! 
    1274                !                                  ! up & down beta terms 
    1275                ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!) 
    1276                IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 
    1277                ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig 
    1278                ENDIF 
    1279                ! 
    1280                IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
    1281                ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig 
    1282                ENDIF 
    1283                ! 
    1284                ! if all the points are outside ice cover 
    1285                IF( zup == -zbig )   zbetup(ji,jj,jl) = 0._wp ! zbig 
    1286                IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0._wp ! zbig             
    1287                ! 
    1288             END DO 
    1289          END DO 
     1149         DO_2D_11_11 
     1150            IF    ( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 
     1151               zbup(ji,jj) = -zbig 
     1152               zbdo(ji,jj) =  zbig 
     1153            ELSEIF( pt(ji,jj,jl) <= 0._wp .AND. pt_ups(ji,jj,jl) > 0._wp ) THEN 
     1154               zbup(ji,jj) = pt_ups(ji,jj,jl) 
     1155               zbdo(ji,jj) = pt_ups(ji,jj,jl) 
     1156            ELSEIF( pt(ji,jj,jl) > 0._wp .AND. pt_ups(ji,jj,jl) <= 0._wp ) THEN 
     1157               zbup(ji,jj) = pt(ji,jj,jl) 
     1158               zbdo(ji,jj) = pt(ji,jj,jl) 
     1159            ELSE 
     1160               zbup(ji,jj) = MAX( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 
     1161               zbdo(ji,jj) = MIN( pt(ji,jj,jl) , pt_ups(ji,jj,jl) ) 
     1162            ENDIF 
     1163         END_2D 
     1164 
     1165         DO_2D_00_00 
     1166            ! 
     1167            zup  = MAX( zbup(ji,jj), zbup(ji-1,jj), zbup(ji+1,jj), zbup(ji,jj-1), zbup(ji,jj+1) )  ! search max/min in neighbourhood 
     1168            zdo  = MIN( zbdo(ji,jj), zbdo(ji-1,jj), zbdo(ji+1,jj), zbdo(ji,jj-1), zbdo(ji,jj+1) ) 
     1169            ! 
     1170            zpos = MAX( 0._wp, pfu_ho(ji-1,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji  ,jj  ,jl) ) &  ! positive/negative part of the flux 
     1171               & + MAX( 0._wp, pfv_ho(ji  ,jj-1,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj  ,jl) ) 
     1172            zneg = MAX( 0._wp, pfu_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfu_ho(ji-1,jj  ,jl) ) & 
     1173               & + MAX( 0._wp, pfv_ho(ji  ,jj  ,jl) ) - MIN( 0._wp, pfv_ho(ji  ,jj-1,jl) ) 
     1174            ! 
     1175            zpos = zpos - (pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1176               &          ) * ( 1. - pamsk ) 
     1177            zneg = zneg + (pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1) ) & 
     1178               &          ) * ( 1. - pamsk ) 
     1179            ! 
     1180            !                                  ! up & down beta terms 
     1181            ! clem: zbetup and zbetdo must be 0 for zpos>1.e-10 & zneg>1.e-10 (do not put 0 instead of 1.e-10 !!!) 
     1182            IF( zpos > epsi10 ) THEN ; zbetup(ji,jj,jl) = MAX( 0._wp, zup - pt_ups(ji,jj,jl) ) / zpos * e1e2t(ji,jj) * z1_dt 
     1183            ELSE                     ; zbetup(ji,jj,jl) = 0._wp ! zbig 
     1184            ENDIF 
     1185            ! 
     1186            IF( zneg > epsi10 ) THEN ; zbetdo(ji,jj,jl) = MAX( 0._wp, pt_ups(ji,jj,jl) - zdo ) / zneg * e1e2t(ji,jj) * z1_dt 
     1187            ELSE                     ; zbetdo(ji,jj,jl) = 0._wp ! zbig 
     1188            ENDIF 
     1189            ! 
     1190            ! if all the points are outside ice cover 
     1191            IF( zup == -zbig )   zbetup(ji,jj,jl) = 0._wp ! zbig 
     1192            IF( zdo ==  zbig )   zbetdo(ji,jj,jl) = 0._wp ! zbig             
     1193            ! 
     1194         END_2D 
    12901195      END DO 
    12911196      CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. )   ! lateral boundary cond. (unchanged sign) 
    … …  
    12951200      ! --------------------------------- 
    12961201      DO jl = 1, jpl 
    1297          DO jj = 1, jpjm1 
    1298             DO ji = 1, fs_jpim1   ! vector opt. 
    1299                zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 
    1300                zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 
    1301                zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) ) 
    1302                ! 
    1303                zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 
    1304                ! 
    1305                pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl) 
    1306                ! 
    1307             END DO 
    1308          END DO 
    1309  
    1310          DO jj = 1, jpjm1 
    1311             DO ji = 1, fs_jpim1   ! vector opt. 
    1312                zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 
    1313                zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 
    1314                zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) ) 
    1315                ! 
    1316                zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 
    1317                ! 
    1318                pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl) 
    1319                ! 
    1320             END DO 
    1321          END DO 
     1202         DO_2D_10_10 
     1203            zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 
     1204            zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 
     1205            zcu = 0.5_wp + SIGN( 0.5_wp , pfu_ho(ji,jj,jl) ) 
     1206            ! 
     1207            zcoef = ( zcu * zau + ( 1._wp - zcu ) * zbu ) 
     1208            ! 
     1209            pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) * zcoef + pfu_ups(ji,jj,jl) 
     1210            ! 
     1211         END_2D 
     1212 
     1213         DO_2D_10_10 
     1214            zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 
     1215            zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 
     1216            zcv = 0.5_wp + SIGN( 0.5_wp , pfv_ho(ji,jj,jl) ) 
     1217            ! 
     1218            zcoef = ( zcv * zav + ( 1._wp - zcv ) * zbv ) 
     1219            ! 
     1220            pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) * zcoef + pfv_ups(ji,jj,jl) 
     1221            ! 
     1222         END_2D 
    13221223 
    13231224      END DO 
    … …  
    13441245      ! 
    13451246      DO jl = 1, jpl 
    1346          DO jj = 2, jpjm1 
    1347             DO ji = fs_2, fs_jpim1   ! vector opt. 
    1348                zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 
    1349             END DO 
    1350          END DO 
     1247         DO_2D_00_00 
     1248            zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 
     1249         END_2D 
    13511250      END DO 
    13521251      CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.)   ! lateral boundary cond. 
    13531252       
    13541253      DO jl = 1, jpl 
    1355          DO jj = 2, jpjm1 
    1356             DO ji = fs_2, fs_jpim1   ! vector opt. 
    1357                uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 
    1358                 
    1359                Rjm = zslpx(ji-1,jj,jl) 
    1360                Rj  = zslpx(ji  ,jj,jl) 
    1361                Rjp = zslpx(ji+1,jj,jl) 
    1362  
    1363                IF( np_limiter == 3 ) THEN 
    1364  
    1365                   IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
    1366                   ELSE                        ;   Rr = Rjp 
     1254         DO_2D_00_00 
     1255            uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 
     1256             
     1257            Rjm = zslpx(ji-1,jj,jl) 
     1258            Rj  = zslpx(ji  ,jj,jl) 
     1259            Rjp = zslpx(ji+1,jj,jl) 
     1260 
     1261            IF( np_limiter == 3 ) THEN 
     1262 
     1263               IF( pu(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     1264               ELSE                        ;   Rr = Rjp 
     1265               ENDIF 
     1266 
     1267               zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)      
     1268               IF( Rj > 0. ) THEN 
     1269                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  & 
     1270                     &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 
     1271               ELSE 
     1272                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  & 
     1273                     &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 
     1274               ENDIF 
     1275               pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 
     1276 
     1277            ELSEIF( np_limiter == 2 ) THEN 
     1278               IF( Rj /= 0. ) THEN 
     1279                  IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     1280                  ELSE                        ;   Cr = Rjp / Rj 
    13671281                  ENDIF 
    1368  
    1369                   zh3 = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)      
    1370                   IF( Rj > 0. ) THEN 
    1371                      zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pu(ji,jj)),  & 
    1372                         &        MIN( 2. * Rr * 0.5 * ABS(pu(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 
    1373                   ELSE 
    1374                      zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pu(ji,jj)),  & 
    1375                         &        MIN(-2. * Rr * 0.5 * ABS(pu(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pu(ji,jj)) ) ) ) ) 
    1376                   ENDIF 
    1377                   pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 
    1378  
    1379                ELSEIF( np_limiter == 2 ) THEN 
    1380                   IF( Rj /= 0. ) THEN 
    1381                      IF( pu(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
    1382                      ELSE                        ;   Cr = Rjp / Rj 
    1383                      ENDIF 
    1384                   ELSE 
    1385                      Cr = 0. 
    1386                   ENDIF 
    1387  
    1388                   ! -- superbee -- 
    1389                   zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
    1390                   ! -- van albada 2 -- 
    1391                   !!zpsi = 2.*Cr / (Cr*Cr+1.) 
    1392                   ! -- sweby (with beta=1) -- 
    1393                   !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
    1394                   ! -- van Leer -- 
    1395                   !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
    1396                   ! -- ospre -- 
    1397                   !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
    1398                   ! -- koren -- 
    1399                   !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
    1400                   ! -- charm -- 
    1401                   !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
    1402                   !ELSE                 ;   zpsi = 0. 
    1403                   !ENDIF 
    1404                   ! -- van albada 1 -- 
    1405                   !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
    1406                   ! -- smart -- 
    1407                   !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
    1408                   ! -- umist -- 
    1409                   !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
    1410  
    1411                   ! high order flux corrected by the limiter 
    1412                   pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 
    1413  
     1282               ELSE 
     1283                  Cr = 0. 
    14141284               ENDIF 
    1415             END DO 
    1416          END DO 
     1285 
     1286               ! -- superbee -- 
     1287               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
     1288               ! -- van albada 2 -- 
     1289               !!zpsi = 2.*Cr / (Cr*Cr+1.) 
     1290               ! -- sweby (with beta=1) -- 
     1291               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
     1292               ! -- van Leer -- 
     1293               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
     1294               ! -- ospre -- 
     1295               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
     1296               ! -- koren -- 
     1297               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
     1298               ! -- charm -- 
     1299               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
     1300               !ELSE                 ;   zpsi = 0. 
     1301               !ENDIF 
     1302               ! -- van albada 1 -- 
     1303               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
     1304               ! -- smart -- 
     1305               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
     1306               ! -- umist -- 
     1307               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
     1308 
     1309               ! high order flux corrected by the limiter 
     1310               pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - ABS( pu(ji,jj) ) * ( (1.-zpsi) + uCFL*zpsi ) * Rj * 0.5 
     1311 
     1312            ENDIF 
     1313         END_2D 
    14171314      END DO 
    14181315      CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.)   ! lateral boundary cond. 
    … …  
    14391336      ! 
    14401337      DO jl = 1, jpl 
    1441          DO jj = 2, jpjm1 
    1442             DO ji = fs_2, fs_jpim1   ! vector opt. 
    1443                zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 
    1444             END DO 
    1445          END DO 
     1338         DO_2D_00_00 
     1339            zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 
     1340         END_2D 
    14461341      END DO 
    14471342      CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.)   ! lateral boundary cond. 
    14481343 
    14491344      DO jl = 1, jpl 
    1450          DO jj = 2, jpjm1 
    1451             DO ji = fs_2, fs_jpim1   ! vector opt. 
    1452                vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 
    1453  
    1454                Rjm = zslpy(ji,jj-1,jl) 
    1455                Rj  = zslpy(ji,jj  ,jl) 
    1456                Rjp = zslpy(ji,jj+1,jl) 
    1457  
    1458                IF( np_limiter == 3 ) THEN 
    1459  
    1460                   IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
    1461                   ELSE                        ;   Rr = Rjp 
     1345         DO_2D_00_00 
     1346            vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 
     1347 
     1348            Rjm = zslpy(ji,jj-1,jl) 
     1349            Rj  = zslpy(ji,jj  ,jl) 
     1350            Rjp = zslpy(ji,jj+1,jl) 
     1351 
     1352            IF( np_limiter == 3 ) THEN 
     1353 
     1354               IF( pv(ji,jj) > 0. ) THEN   ;   Rr = Rjm 
     1355               ELSE                        ;   Rr = Rjp 
     1356               ENDIF 
     1357 
     1358               zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)      
     1359               IF( Rj > 0. ) THEN 
     1360                  zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  & 
     1361                     &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 
     1362               ELSE 
     1363                  zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  & 
     1364                     &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 
     1365               ENDIF 
     1366               pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 
     1367 
     1368            ELSEIF( np_limiter == 2 ) THEN 
     1369 
     1370               IF( Rj /= 0. ) THEN 
     1371                  IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
     1372                  ELSE                        ;   Cr = Rjp / Rj 
    14621373                  ENDIF 
    1463  
    1464                   zh3 = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)      
    1465                   IF( Rj > 0. ) THEN 
    1466                      zlimiter =  MAX( 0., MIN( zh3, MAX(-Rr * 0.5 * ABS(pv(ji,jj)),  & 
    1467                         &        MIN( 2. * Rr * 0.5 * ABS(pv(ji,jj)),  zh3,  1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 
    1468                   ELSE 
    1469                      zlimiter = -MAX( 0., MIN(-zh3, MAX( Rr * 0.5 * ABS(pv(ji,jj)),  & 
    1470                         &        MIN(-2. * Rr * 0.5 * ABS(pv(ji,jj)), -zh3, -1.5 * Rj * 0.5 * ABS(pv(ji,jj)) ) ) ) ) 
    1471                   ENDIF 
    1472                   pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 
    1473  
    1474                ELSEIF( np_limiter == 2 ) THEN 
    1475  
    1476                   IF( Rj /= 0. ) THEN 
    1477                      IF( pv(ji,jj) > 0. ) THEN   ;   Cr = Rjm / Rj 
    1478                      ELSE                        ;   Cr = Rjp / Rj 
    1479                      ENDIF 
    1480                   ELSE 
    1481                      Cr = 0. 
    1482                   ENDIF 
    1483  
    1484                   ! -- superbee -- 
    1485                   zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
    1486                   ! -- van albada 2 -- 
    1487                   !!zpsi = 2.*Cr / (Cr*Cr+1.) 
    1488                   ! -- sweby (with beta=1) -- 
    1489                   !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
    1490                   ! -- van Leer -- 
    1491                   !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
    1492                   ! -- ospre -- 
    1493                   !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
    1494                   ! -- koren -- 
    1495                   !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
    1496                   ! -- charm -- 
    1497                   !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
    1498                   !ELSE                 ;   zpsi = 0. 
    1499                   !ENDIF 
    1500                   ! -- van albada 1 -- 
    1501                   !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
    1502                   ! -- smart -- 
    1503                   !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
    1504                   ! -- umist -- 
    1505                   !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
    1506  
    1507                   ! high order flux corrected by the limiter 
    1508                   pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 
    1509  
     1374               ELSE 
     1375                  Cr = 0. 
    15101376               ENDIF 
    1511             END DO 
    1512          END DO 
     1377 
     1378               ! -- superbee -- 
     1379               zpsi = MAX( 0., MAX( MIN(1.,2.*Cr), MIN(2.,Cr) ) ) 
     1380               ! -- van albada 2 -- 
     1381               !!zpsi = 2.*Cr / (Cr*Cr+1.) 
     1382               ! -- sweby (with beta=1) -- 
     1383               !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) 
     1384               ! -- van Leer -- 
     1385               !!zpsi = ( Cr + ABS(Cr) ) / ( 1. + ABS(Cr) ) 
     1386               ! -- ospre -- 
     1387               !!zpsi = 1.5 * ( Cr*Cr + Cr ) / ( Cr*Cr + Cr + 1. ) 
     1388               ! -- koren -- 
     1389               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( (1.+2*Cr)/3., 2. ) ) ) 
     1390               ! -- charm -- 
     1391               !IF( Cr > 0. ) THEN   ;   zpsi = Cr * (3.*Cr + 1.) / ( (Cr + 1.) * (Cr + 1.) ) 
     1392               !ELSE                 ;   zpsi = 0. 
     1393               !ENDIF 
     1394               ! -- van albada 1 -- 
     1395               !!zpsi = (Cr*Cr + Cr) / (Cr*Cr +1) 
     1396               ! -- smart -- 
     1397               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, 4. ) ) ) 
     1398               ! -- umist -- 
     1399               !!zpsi = MAX( 0., MIN( 2.*Cr, MIN( 0.25+0.75*Cr, MIN(0.75+0.25*Cr, 2. ) ) ) ) 
     1400 
     1401               ! high order flux corrected by the limiter 
     1402               pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - ABS( pv(ji,jj) ) * ( (1.-zpsi) + vCFL*zpsi ) * Rj * 0.5 
     1403 
     1404            ENDIF 
     1405         END_2D 
    15131406      END DO 
    15141407      CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.)   ! lateral boundary cond. 
    … …  
    15441437      DO jl = 1, jpl 
    15451438 
    1546          DO jj = 1, jpj 
    1547             DO ji = 1, jpi 
    1548                IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     1439         DO_2D_11_11 
     1440            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     1441               ! 
     1442               !                               ! -- check h_ip -- ! 
     1443               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     1444               IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1445                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
     1446                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     1447                     pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     1448                  ENDIF 
     1449               ENDIF 
     1450               ! 
     1451               !                               ! -- check h_i -- ! 
     1452               ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     1453               zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
     1454               IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1455                  pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     1456               ENDIF 
     1457               ! 
     1458               !                               ! -- check h_s -- ! 
     1459               ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     1460               zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
     1461               IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1462                  zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
    15491463                  ! 
    1550                   !                               ! -- check h_ip -- ! 
    1551                   ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    1552                   IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    1553                      zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    1554                      IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
    1555                         pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
    1556                      ENDIF 
    1557                   ENDIF 
     1464                  wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 
     1465                  hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
    15581466                  ! 
    1559                   !                               ! -- check h_i -- ! 
    1560                   ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
    1561                   zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 
    1562                   IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
    1563                      pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
    1564                   ENDIF 
    1565                   ! 
    1566                   !                               ! -- check h_s -- ! 
    1567                   ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
    1568                   zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 
    1569                   IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
    1570                      zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 
    1571                      ! 
    1572                      wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 
    1573                      hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
    1574                      ! 
    1575                      pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
    1576                      pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
    1577                   ENDIF            
    1578                   !                   
    1579                ENDIF 
    1580             END DO 
    1581          END DO 
     1467                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     1468                  pv_s(ji,jj,jl)          = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     1469               ENDIF            
     1470               !                   
     1471            ENDIF 
     1472         END_2D 
    15821473      END DO  
    15831474      ! 
    … …  
    16121503      ! -- check snow load -- ! 
    16131504      DO jl = 1, jpl 
    1614          DO jj = 1, jpj 
    1615             DO ji = 1, jpi 
    1616                IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
    1617                   ! 
    1618                   zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
    1619                   ! 
    1620                   IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface 
    1621                      ! put snow excess in the ocean 
    1622                      zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
    1623                      wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 
    1624                      hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
    1625                      ! correct snow volume and heat content 
    1626                      pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
    1627                      pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
    1628                   ENDIF 
    1629                   ! 
     1505         DO_2D_11_11 
     1506            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     1507               ! 
     1508               zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
     1509               ! 
     1510               IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface 
     1511                  ! put snow excess in the ocean 
     1512                  zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 
     1513                  wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 
     1514                  hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1515                  ! correct snow volume and heat content 
     1516                  pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 
     1517                  pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess 
    16301518               ENDIF 
    1631             END DO 
    1632          END DO 
     1519               ! 
     1520            ENDIF 
     1521         END_2D 
    16331522      END DO 
    16341523      ! 
Note: See TracChangeset for help on using the changeset viewer.