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

    r12236 r12340  
    4949   !! * Substitutions 
    5050#  include "vectopt_loop_substitute.h90" 
     51#  include "do_loop_substitute.h90" 
    5152   !!---------------------------------------------------------------------- 
    5253   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    … …  
    180181      !------------------------------------------------------------------------------! 
    181182      ! ocean/land mask 
    182       DO jj = 1, jpjm1 
    183          DO ji = 1, jpim1      ! NO vector opt. 
    184             zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 
    185          END DO 
    186       END DO 
     183      DO_2D_10_10 
     184         zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 
     185      END_2D 
    187186      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp ) 
    188187 
    189188      ! Lateral boundary conditions on velocity (modify zfmask) 
    190189      zwf(:,:) = zfmask(:,:) 
    191       DO jj = 2, jpjm1 
    192          DO ji = fs_2, fs_jpim1   ! vector opt. 
    193             IF( zfmask(ji,jj) == 0._wp ) THEN 
    194                zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 
    195             ENDIF 
    196          END DO 
    197       END DO 
     190      DO_2D_00_00 
     191         IF( zfmask(ji,jj) == 0._wp ) THEN 
     192            zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 
     193         ENDIF 
     194      END_2D 
    198195      DO jj = 2, jpjm1 
    199196         IF( zfmask(1,jj) == 0._wp ) THEN 
    … …  
    257254      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
    258255 
    259       DO jj = 2, jpjm1 
    260          DO ji = fs_2, fs_jpim1 
    261  
    262             ! ice fraction at U-V points 
    263             zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    264             zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    265  
    266             ! Ice/snow mass at U-V points 
    267             zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) ) 
    268             zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) ) 
    269             zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) ) 
    270             zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    271             zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    272  
    273             ! Ocean currents at U-V points 
    274             v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 
    275             u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 
    276  
    277             ! Coriolis at T points (m*f) 
    278             zmf(ji,jj)      = zm1 * ff_t(ji,jj) 
    279  
    280             ! dt/m at T points (for alpha and beta coefficients) 
    281             zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin ) 
    282              
    283             ! m/dt 
    284             zmU_t(ji,jj)    = zmassU * z1_dtevp 
    285             zmV_t(ji,jj)    = zmassV * z1_dtevp 
    286              
    287             ! Drag ice-atm. 
    288             ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
    289             ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
    290  
    291             ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
    292             zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 
    293             zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 
    294  
    295             ! masks 
    296             zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
    297             zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
    298  
    299             ! switches 
    300             IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
    301             ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
    302             IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
    303             ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF 
    304  
    305          END DO 
    306       END DO 
     256      DO_2D_00_00 
     257 
     258         ! ice fraction at U-V points 
     259         zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
     260         zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     261 
     262         ! Ice/snow mass at U-V points 
     263         zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) ) 
     264         zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) ) 
     265         zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) ) 
     266         zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
     267         zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     268 
     269         ! Ocean currents at U-V points 
     270         v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 
     271         u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 
     272 
     273         ! Coriolis at T points (m*f) 
     274         zmf(ji,jj)      = zm1 * ff_t(ji,jj) 
     275 
     276         ! dt/m at T points (for alpha and beta coefficients) 
     277         zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin ) 
     278          
     279         ! m/dt 
     280         zmU_t(ji,jj)    = zmassU * z1_dtevp 
     281         zmV_t(ji,jj)    = zmassV * z1_dtevp 
     282          
     283         ! Drag ice-atm. 
     284         ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     285         ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
     286 
     287         ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
     288         zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj) 
     289         zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj) 
     290 
     291         ! masks 
     292         zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
     293         zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
     294 
     295         ! switches 
     296         IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
     297         ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
     298         IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
     299         ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF 
     300 
     301      END_2D 
    307302      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. ) 
    308303      ! 
    … …  
    310305      ! 
    311306      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016 
    312          DO jj = 2, jpjm1 
    313             DO ji = fs_2, fs_jpim1 
    314                ! ice thickness at U-V points 
    315                zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    316                zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    317                ! ice-bottom stress at U points 
    318                zvCr = zaU(ji,jj) * rn_depfra * hu(ji,jj,Kmm) 
    319                ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    320                ! ice-bottom stress at V points 
    321                zvCr = zaV(ji,jj) * rn_depfra * hv(ji,jj,Kmm) 
    322                ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    323                ! ice_bottom stress at T points 
    324                zvCr = at_i(ji,jj) * rn_depfra * ht(ji,jj) 
    325                tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    326             END DO 
    327          END DO 
     307         DO_2D_00_00 
     308            ! ice thickness at U-V points 
     309            zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
     310            zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     311            ! ice-bottom stress at U points 
     312            zvCr = zaU(ji,jj) * rn_depfra * hu(ji,jj,Kmm) 
     313            ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
     314            ! ice-bottom stress at V points 
     315            zvCr = zaV(ji,jj) * rn_depfra * hv(ji,jj,Kmm) 
     316            ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
     317            ! ice_bottom stress at T points 
     318            zvCr = at_i(ji,jj) * rn_depfra * ht(ji,jj) 
     319            tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
     320         END_2D 
    328321         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 
    329322         ! 
    330323      ELSE                               !-- no landfast 
    331          DO jj = 2, jpjm1 
    332             DO ji = fs_2, fs_jpim1 
    333                ztaux_base(ji,jj) = 0._wp 
    334                ztauy_base(ji,jj) = 0._wp 
    335             END DO 
    336          END DO 
     324         DO_2D_00_00 
     325            ztaux_base(ji,jj) = 0._wp 
     326            ztauy_base(ji,jj) = 0._wp 
     327         END_2D 
    337328      ENDIF 
    338329 
    … …  
    354345 
    355346         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
    356          DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
    357             DO ji = 1, jpim1 
    358  
    359                ! shear at F points 
    360                zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    361                   &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    362                   &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
    363  
    364             END DO 
    365          END DO 
     347         DO_2D_10_10 
     348 
     349            ! shear at F points 
     350            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     351               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     352               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
     353 
     354         END_2D 
    366355         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. ) 
    367356 
    368          DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
    369             DO ji = 2, jpi ! no vector loop 
    370  
    371                ! shear**2 at T points (doc eq. A16) 
    372                zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
    373                   &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
    374                   &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
    375                
    376                ! divergence at T points 
    377                zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    378                   &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    379                   &    ) * r1_e1e2t(ji,jj) 
    380                zdiv2 = zdiv * zdiv 
    381                 
    382                ! tension at T points 
    383                zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    384                   &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    385                   &   ) * r1_e1e2t(ji,jj) 
    386                zdt2 = zdt * zdt 
    387                 
    388                ! delta at T points 
    389                zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
    390  
    391                ! P/delta at T points 
    392                zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
    393  
    394                ! alpha & beta for aEVP 
    395                !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
    396                !   alpha = beta = sqrt(4*gamma) 
    397                IF( ln_aEVP ) THEN 
    398                   zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
    399                   z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    400                   zalph2   = zalph1 
    401                   z1_alph2 = z1_alph1 
    402                ENDIF 
    403                 
    404                ! stress at T points (zkt/=0 if landfast) 
    405                zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
    406                zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
    407               
    408             END DO 
    409          END DO 
     357         DO_2D_01_01 
     358 
     359            ! shear**2 at T points (doc eq. A16) 
     360            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     361               &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
     362               &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
     363            
     364            ! divergence at T points 
     365            zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     366               &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     367               &    ) * r1_e1e2t(ji,jj) 
     368            zdiv2 = zdiv * zdiv 
     369             
     370            ! tension at T points 
     371            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     372               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     373               &   ) * r1_e1e2t(ji,jj) 
     374            zdt2 = zdt * zdt 
     375             
     376            ! delta at T points 
     377            zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     378 
     379            ! P/delta at T points 
     380            zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
     381 
     382            ! alpha & beta for aEVP 
     383            !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
     384            !   alpha = beta = sqrt(4*gamma) 
     385            IF( ln_aEVP ) THEN 
     386               zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     387               z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     388               zalph2   = zalph1 
     389               z1_alph2 = z1_alph1 
     390            ENDIF 
     391             
     392            ! stress at T points (zkt/=0 if landfast) 
     393            zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
     394            zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
     395           
     396         END_2D 
    410397         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. ) 
    411398 
    412          DO jj = 1, jpjm1 
    413             DO ji = 1, jpim1 
    414  
    415                ! alpha & beta for aEVP 
    416                IF( ln_aEVP ) THEN 
    417                   zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
    418                   z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
    419                   zbeta(ji,jj) = zalph2 
    420                ENDIF 
    421                 
    422                ! P/delta at F points 
    423                zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
    424                 
    425                ! stress at F points (zkt/=0 if landfast) 
    426                zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 
    427  
    428             END DO 
    429          END DO 
     399         DO_2D_10_10 
     400 
     401            ! alpha & beta for aEVP 
     402            IF( ln_aEVP ) THEN 
     403               zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     404               z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     405               zbeta(ji,jj) = zalph2 
     406            ENDIF 
     407             
     408            ! P/delta at F points 
     409            zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
     410             
     411            ! stress at F points (zkt/=0 if landfast) 
     412            zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 
     413 
     414         END_2D 
    430415 
    431416         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
    432          DO jj = 2, jpjm1 
    433             DO ji = fs_2, fs_jpim1                
    434                !                   !--- U points 
    435                zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
    436                   &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
    437                   &                    ) * r1_e2u(ji,jj)                                                                      & 
    438                   &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
    439                   &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
    440                   &                  ) * r1_e1e2u(ji,jj) 
    441                ! 
    442                !                !--- V points 
    443                zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
    444                   &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
    445                   &                    ) * r1_e1v(ji,jj)                                                                      & 
    446                   &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
    447                   &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
    448                   &                  ) * r1_e1e2v(ji,jj) 
    449                ! 
    450                !                !--- ice currents at U-V point 
    451                v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) 
    452                u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) 
    453                ! 
    454             END DO 
    455          END DO 
     417         DO_2D_00_00 
     418            !                   !--- U points 
     419            zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
     420               &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     421               &                    ) * r1_e2u(ji,jj)                                                                      & 
     422               &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
     423               &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
     424               &                  ) * r1_e1e2u(ji,jj) 
     425            ! 
     426            !                !--- V points 
     427            zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
     428               &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     429               &                    ) * r1_e1v(ji,jj)                                                                      & 
     430               &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
     431               &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
     432               &                  ) * r1_e1e2v(ji,jj) 
     433            ! 
     434            !                !--- ice currents at U-V point 
     435            v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) 
     436            u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) 
     437            ! 
     438         END_2D 
    456439         ! 
    457440         ! --- Computation of ice velocity --- ! 
    … …  
    460443         IF( MOD(jter,2) == 0 ) THEN ! even iterations 
    461444            ! 
    462             DO jj = 2, jpjm1 
    463                DO ji = fs_2, fs_jpim1 
    464                   !                 !--- tau_io/(v_oce - v_ice) 
    465                   zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    466                      &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    467                   !                 !--- Ocean-to-Ice stress 
    468                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
    469                   ! 
    470                   !                 !--- tau_bottom/v_ice 
    471                   zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    472                   zTauB = ztauy_base(ji,jj) / zvel 
    473                   !                 !--- OceanBottom-to-Ice stress 
    474                   ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    475                   ! 
    476                   !                 !--- Coriolis at V-points (energy conserving formulation) 
    477                   zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    478                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    479                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    480                   ! 
    481                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    482                   zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    483                   ! 
    484                   !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
    485                   !                                         1 = sliding friction : TauB < RHS 
    486                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    487                   ! 
    488                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    489                      v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    490                         &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    491                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    492                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    493                         &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    494                         &           )   * zmsk00y(ji,jj) 
    495                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    496                      v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
    497                         &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    498                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    499                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    500                         &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    501                         &            )   * zmsk00y(ji,jj) 
    502                   ENDIF 
    503                END DO 
    504             END DO 
     445            DO_2D_00_00 
     446               !                 !--- tau_io/(v_oce - v_ice) 
     447               zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     448                  &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     449               !                 !--- Ocean-to-Ice stress 
     450               ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     451               ! 
     452               !                 !--- tau_bottom/v_ice 
     453               zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
     454               zTauB = ztauy_base(ji,jj) / zvel 
     455               !                 !--- OceanBottom-to-Ice stress 
     456               ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
     457               ! 
     458               !                 !--- Coriolis at V-points (energy conserving formulation) 
     459               zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     460                  &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     461                  &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     462               ! 
     463               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     464               zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     465               ! 
     466               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     467               !                                         1 = sliding friction : TauB < RHS 
     468               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
     469               ! 
     470               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     471                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
     472                     &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     473                     &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     474                     &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     475                     &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     476                     &           )   * zmsk00y(ji,jj) 
     477               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
     478                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
     479                     &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     480                     &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     481                     &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     482                     &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     483                     &            )   * zmsk00y(ji,jj) 
     484               ENDIF 
     485            END_2D 
    505486            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. ) 
    506487            ! 
    … …  
    511492            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' ) 
    512493            ! 
    513             DO jj = 2, jpjm1 
    514                DO ji = fs_2, fs_jpim1           
    515                   !                 !--- tau_io/(u_oce - u_ice) 
    516                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    517                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    518                   !                 !--- Ocean-to-Ice stress 
    519                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    520                   ! 
    521                   !                 !--- tau_bottom/u_ice 
    522                   zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    523                   zTauB = ztaux_base(ji,jj) / zvel 
    524                   !                 !--- OceanBottom-to-Ice stress 
    525                   ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
    526                   ! 
    527                   !                 !--- Coriolis at U-points (energy conserving formulation) 
    528                   zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    529                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    530                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    531                   ! 
    532                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    533                   zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    534                   ! 
    535                   !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
    536                   !                                         1 = sliding friction : TauB < RHS 
    537                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    538                   ! 
    539                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    540                      u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    541                         &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    542                         &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    543                         &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    544                         &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    545                         &           )   * zmsk00x(ji,jj) 
    546                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    547                      u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
    548                         &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    549                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    550                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    551                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    552                         &            )   * zmsk00x(ji,jj) 
    553                   ENDIF 
    554                END DO 
    555             END DO 
     494            DO_2D_00_00 
     495               !                 !--- tau_io/(u_oce - u_ice) 
     496               zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     497                  &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     498               !                 !--- Ocean-to-Ice stress 
     499               ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     500               ! 
     501               !                 !--- tau_bottom/u_ice 
     502               zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
     503               zTauB = ztaux_base(ji,jj) / zvel 
     504               !                 !--- OceanBottom-to-Ice stress 
     505               ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
     506               ! 
     507               !                 !--- Coriolis at U-points (energy conserving formulation) 
     508               zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
     509                  &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     510                  &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     511               ! 
     512               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     513               zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     514               ! 
     515               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     516               !                                         1 = sliding friction : TauB < RHS 
     517               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
     518               ! 
     519               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     520                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
     521                     &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     522                     &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     523                     &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     524                     &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     525                     &           )   * zmsk00x(ji,jj) 
     526               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
     527                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
     528                     &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     529                     &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     530                     &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     531                     &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     532                     &            )   * zmsk00x(ji,jj) 
     533               ENDIF 
     534            END_2D 
    556535            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. ) 
    557536            ! 
    … …  
    564543         ELSE ! odd iterations 
    565544            ! 
    566             DO jj = 2, jpjm1 
    567                DO ji = fs_2, fs_jpim1 
    568                   !                 !--- tau_io/(u_oce - u_ice) 
    569                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    570                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    571                   !                 !--- Ocean-to-Ice stress 
    572                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    573                   ! 
    574                   !                 !--- tau_bottom/u_ice 
    575                   zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    576                   zTauB = ztaux_base(ji,jj) / zvel 
    577                   !                 !--- OceanBottom-to-Ice stress 
    578                   ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
    579                   ! 
    580                   !                 !--- Coriolis at U-points (energy conserving formulation) 
    581                   zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    582                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    583                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    584                   ! 
    585                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    586                   zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    587                   ! 
    588                   !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
    589                   !                                         1 = sliding friction : TauB < RHS 
    590                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    591                   ! 
    592                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    593                      u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    594                         &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    595                         &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    596                         &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    597                         &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    598                         &           )   * zmsk00x(ji,jj) 
    599                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    600                      u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
    601                         &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    602                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    603                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    604                         &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    605                         &            )   * zmsk00x(ji,jj) 
    606                   ENDIF 
    607                END DO 
    608             END DO 
     545            DO_2D_00_00 
     546               !                 !--- tau_io/(u_oce - u_ice) 
     547               zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     548                  &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     549               !                 !--- Ocean-to-Ice stress 
     550               ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     551               ! 
     552               !                 !--- tau_bottom/u_ice 
     553               zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
     554               zTauB = ztaux_base(ji,jj) / zvel 
     555               !                 !--- OceanBottom-to-Ice stress 
     556               ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
     557               ! 
     558               !                 !--- Coriolis at U-points (energy conserving formulation) 
     559               zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
     560                  &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     561                  &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     562               ! 
     563               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     564               zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     565               ! 
     566               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     567               !                                         1 = sliding friction : TauB < RHS 
     568               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
     569               ! 
     570               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     571                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
     572                     &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     573                     &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     574                     &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     575                     &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     576                     &           )   * zmsk00x(ji,jj) 
     577               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
     578                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
     579                     &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     580                     &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     581                     &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     582                     &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     583                     &            )   * zmsk00x(ji,jj) 
     584               ENDIF 
     585            END_2D 
    609586            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. ) 
    610587            ! 
    … …  
    615592            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' ) 
    616593            ! 
    617             DO jj = 2, jpjm1 
    618                DO ji = fs_2, fs_jpim1 
    619                   !                 !--- tau_io/(v_oce - v_ice) 
    620                   zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    621                      &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    622                   !                 !--- Ocean-to-Ice stress 
    623                   ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
    624                   ! 
    625                   !                 !--- tau_bottom/v_ice 
    626                   zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    627                   zTauB = ztauy_base(ji,jj) / zvel 
    628                   !                 !--- OceanBottom-to-Ice stress 
    629                   ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    630                   ! 
    631                   !                 !--- Coriolis at v-points (energy conserving formulation) 
    632                   zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    633                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    634                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    635                   ! 
    636                   !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    637                   zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    638                   ! 
    639                   !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
    640                   !                                         1 = sliding friction : TauB < RHS 
    641                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    642                   ! 
    643                   IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    644                      v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    645                         &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    646                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    647                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    648                         &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    649                         &           )   * zmsk00y(ji,jj) 
    650                   ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    651                      v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
    652                         &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    653                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    654                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    655                         &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    656                         &            )   * zmsk00y(ji,jj) 
    657                   ENDIF 
    658                END DO 
    659             END DO 
     594            DO_2D_00_00 
     595               !                 !--- tau_io/(v_oce - v_ice) 
     596               zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     597                  &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     598               !                 !--- Ocean-to-Ice stress 
     599               ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     600               ! 
     601               !                 !--- tau_bottom/v_ice 
     602               zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
     603               zTauB = ztauy_base(ji,jj) / zvel 
     604               !                 !--- OceanBottom-to-Ice stress 
     605               ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
     606               ! 
     607               !                 !--- Coriolis at v-points (energy conserving formulation) 
     608               zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     609                  &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     610                  &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     611               ! 
     612               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     613               zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     614               ! 
     615               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     616               !                                         1 = sliding friction : TauB < RHS 
     617               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
     618               ! 
     619               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     620                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
     621                     &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     622                     &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     623                     &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     624                     &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     625                     &           )   * zmsk00y(ji,jj) 
     626               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
     627                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
     628                     &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     629                     &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     630                     &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     631                     &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     632                     &            )   * zmsk00y(ji,jj) 
     633               ENDIF 
     634            END_2D 
    660635            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. ) 
    661636            ! 
    … …  
    683658      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)  
    684659      !------------------------------------------------------------------------------! 
    685       DO jj = 1, jpjm1 
    686          DO ji = 1, jpim1 
    687  
    688             ! shear at F points 
    689             zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    690                &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    691                &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
    692  
    693          END DO 
    694       END DO            
     660      DO_2D_10_10 
     661 
     662         ! shear at F points 
     663         zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     664            &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     665            &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
     666 
     667      END_2D 
    695668       
    696       DO jj = 2, jpjm1 
    697          DO ji = 2, jpim1 ! no vector loop 
    698              
    699             ! tension**2 at T points 
    700             zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    701                &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    702                &   ) * r1_e1e2t(ji,jj) 
    703             zdt2 = zdt * zdt 
    704              
    705             ! shear**2 at T points (doc eq. A16) 
    706             zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
    707                &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
    708                &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
    709              
    710             ! shear at T points 
    711             pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
    712  
    713             ! divergence at T points 
    714             pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    715                &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    716                &             ) * r1_e1e2t(ji,jj) 
    717              
    718             ! delta at T points 
    719             zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    720             rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    721             pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
    722  
    723          END DO 
    724       END DO 
     669      DO_2D_00_00 
     670          
     671         ! tension**2 at T points 
     672         zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     673            &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     674            &   ) * r1_e1e2t(ji,jj) 
     675         zdt2 = zdt * zdt 
     676          
     677         ! shear**2 at T points (doc eq. A16) 
     678         zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     679            &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
     680            &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
     681          
     682         ! shear at T points 
     683         pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     684 
     685         ! divergence at T points 
     686         pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     687            &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     688            &             ) * r1_e1e2t(ji,jj) 
     689          
     690         ! delta at T points 
     691         zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
     692         rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
     693         pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     694 
     695      END_2D 
    725696      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
    726697       
    … …  
    735706      ! 5) diagnostics 
    736707      !------------------------------------------------------------------------------! 
    737       DO jj = 1, jpj 
    738          DO ji = 1, jpi 
    739             zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    740          END DO 
    741       END DO 
     708      DO_2D_11_11 
     709         zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
     710      END_2D 
    742711 
    743712      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
    … …  
    766735         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
    767736         !          
    768          DO jj = 2, jpjm1 
    769             DO ji = 2, jpim1 
    770                zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    771                   &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
    772                   &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
    773  
    774                zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    775  
    776                zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
     737         DO_2D_00_00 
     738            zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
     739               &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
     740               &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
     741 
     742            zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
     743 
     744            zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    777745 
    778746!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    … …  
    780748!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
    781749!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    782                zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    783                zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
    784                zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    785             END DO 
    786          END DO 
     750            zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
     751            zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
     752            zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
     753         END_2D 
    787754         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
    788755         ! 
    … …  
    819786            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 
    820787         ! 
    821          DO jj = 2, jpjm1 
    822             DO ji = 2, jpim1 
    823                ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    824                zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 
    825                zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 
    826  
    827                zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
    828                zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
    829  
    830                zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
    831                zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
    832  
    833                zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component 
    834                zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   '' 
    835  
    836             END DO 
    837          END DO 
     788         DO_2D_00_00 
     789            ! 2D ice mass, snow mass, area transport arrays (X, Y) 
     790            zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 
     791            zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 
     792 
     793            zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
     794            zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
     795 
     796            zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
     797            zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
     798 
     799            zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component 
     800            zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   '' 
     801 
     802         END_2D 
    838803 
    839804         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 
Note: See TracChangeset for help on using the changeset viewer.