New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 12377 for NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 – NEMO

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

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

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

svn merge --ignore-ancestry \

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

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

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

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

Location:
NEMO/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk

    • Property svn:externals

      •  

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

  • NEMO/trunk/src/ICE/icedyn_rhg_evp.F90

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