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 7646 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 – NEMO

Ignore:
Timestamp:
2017-02-06T10:25:03+01:00 (7 years ago)
Author:
timgraham
Message:

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge ​svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r6416 r7646  
    99   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas 
    1010   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation  
    11    !!            3.5  !  2012-08  (R. Benshila)  AGRIF  
     11   !!            3.5  !  2012-08  (R. Benshila)  AGRIF 
     12   !!            3.6  !  2016-06  (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013) 
    1213   !!---------------------------------------------------------------------- 
    13 #if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp ) 
     14#if defined key_lim3 
    1415   !!---------------------------------------------------------------------- 
    15    !!   'key_lim3'               OR                     LIM-3 sea-ice model 
    16    !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model 
     16   !!   'key_lim3'                                      LIM-3 sea-ice model 
    1717   !!---------------------------------------------------------------------- 
    1818   !!   lim_rhg       : computes ice velocities 
    … …  
    2424   USE sbc_oce        ! Surface boundary condition: ocean fields 
    2525   USE sbc_ice        ! Surface boundary condition: ice fields 
    26 #if defined key_lim3 
    27    USE ice            ! LIM-3: ice variables 
    28    USE dom_ice        ! LIM-3: ice domain 
    29    USE limitd_me      ! LIM-3:  
    30 #else 
    31    USE ice_2          ! LIM-2: ice variables 
    32    USE dom_ice_2      ! LIM-2: ice domain 
    33 #endif 
     26   USE ice            ! ice variables 
     27   USE limitd_me      ! ice strength 
    3428   USE lbclnk         ! Lateral Boundary Condition / MPP link 
    3529   USE lib_mpp        ! MPP library 
    … …  
    3832   USE prtctl         ! Print control 
    3933   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    40 #if defined key_agrif && defined key_lim2 
    41    USE agrif_lim2_interp 
     34#if defined key_agrif 
     35   USE agrif_lim3_interp 
    4236#endif 
    43 #if defined key_bdy 
    44    USE bdyice_lim 
    45 #endif 
     37   USE bdy_oce   , ONLY: ln_bdy  
     38   USE bdyice_lim  
    4639 
    4740   IMPLICIT NONE 
    4841   PRIVATE 
    4942 
    50    PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2) 
     43   PUBLIC   lim_rhg        ! routine called by lim_dyn 
    5144 
    5245   !! * Substitutions 
    … …  
    5952CONTAINS 
    6053 
    61    SUBROUTINE lim_rhg( k_j1, k_jpj ) 
     54   SUBROUTINE lim_rhg 
    6255      !!------------------------------------------------------------------- 
    6356      !!                 ***  SUBROUTINE lim_rhg  *** 
    … …  
    10699      !!                 e.g. in the Canadian Archipelago 
    107100      !! 
     101      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013 
     102      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters) 
     103      !!              but this solution appears very unstable (see Kimmritz et al 2016) 
     104      !! 
    108105      !! References : Hunke and Dukowicz, JPO97 
    109106      !!              Bouillon et al., Ocean Modelling 2009 
     107      !!              Bouillon et al., Ocean Modelling 2013 
    110108      !!------------------------------------------------------------------- 
    111       INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation 
    112       INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation 
    113       !! 
    114       INTEGER ::   ji, jj   ! dummy loop indices 
    115       INTEGER ::   jter     ! local integers 
     109      INTEGER ::   ji, jj       ! dummy loop indices 
     110      INTEGER ::   jter         ! local integers 
    116111      CHARACTER (len=50) ::   charout 
    117       REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         ! 
    118       REAL(wp) ::   za, zstms          ! local scalars 
    119       REAL(wp) ::   zc1, zc2, zc3      ! ice mass 
    120  
    121       REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling 
    122       REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 
    123       REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars 
    124       REAL(wp) ::   zu_ice2, zv_ice1              ! 
    125       REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre 
    126       REAL(wp) ::   zdst                          ! shear at the center of the grid point 
    127       REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface 
    128       REAL(wp) ::   sigma1, sigma2                ! internal ice stress 
    129  
    130       REAL(wp) ::   zresm         ! Maximal error on ice velocity 
    131       REAL(wp) ::   zintb, zintn  ! dummy argument 
    132  
    133       REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength 
    134       REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc) 
    135       REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points 
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points 
    137       REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points 
    138       REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays 
    139       REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                            
    140       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points 
    141       REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point 
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses 
    143       REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points 
     112 
     113      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio 
     114      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling 
     115      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity 
     116      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013 
     117      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass 
     118      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars 
     119      REAL(wp) ::   zTauO, zTauB, zTauE, zCor, zvel                          ! temporary scalars 
     120 
     121      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress 
     122      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity 
     123      REAL(wp) ::   zintb, zintn                                             ! dummy argument 
    144124       
    145       REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells 
    146       REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells 
    147       REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2  
    148       REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12 
    149       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity 
    150       REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope: 
    151                                                               !   ocean surface (ssh_m) if ice is not embedded 
    152                                                               !   ice top surface if ice is embedded    
    153  
    154       REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter 
    155       REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 
     125      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors 
     126      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                         ! P/delta at T points 
     127      ! 
     128      REAL(wp), POINTER, DIMENSION(:,:) ::   zaU   , zaV                     ! ice fraction on U/V points 
     129      REAL(wp), POINTER, DIMENSION(:,:) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points 
     130      REAL(wp), POINTER, DIMENSION(:,:) ::   zmf                             ! coriolis parameter at T points 
     131      REAL(wp), POINTER, DIMENSION(:,:) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
     132      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points 
     133      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
     134      REAL(wp), POINTER, DIMENSION(:,:) ::   zfU   , zfV                     ! internal stresses 
     135       
     136      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                             ! shear 
     137      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12                  ! stress tensor components 
     138      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr           ! check convergence 
     139      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice                           ! array used for the calculation of ice surface slope: 
     140                                                                             !   ocean surface (ssh_m) if ice is not embedded 
     141                                                                             !   ice top surface if ice is embedded    
     142      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays 
     143      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence 
     144      REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice 
     145 
     146      REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
     147      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity 
     148      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip) 
    156149      !!------------------------------------------------------------------- 
    157150 
    158       CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    159       CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    160       CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    161       CALL wrk_alloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 ) 
    162  
    163 #if  defined key_lim2 && ! defined key_lim2_vp 
    164 # if defined key_agrif 
    165       USE ice_2, vt_s => hsnm 
    166       USE ice_2, vt_i => hicm 
    167 # else 
    168       vt_s => hsnm 
    169       vt_i => hicm 
    170 # endif 
    171       at_i(:,:) = 1. - frld(:,:) 
    172 #endif 
    173 #if defined key_agrif && defined key_lim2  
    174       CALL agrif_rhg_lim2_load      ! First interpolation of coarse values 
     151      CALL wrk_alloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt ) 
     152      CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 
     153      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
     154      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     155      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
     156 
     157#if defined key_agrif  
     158      CALL agrif_interp_lim3( 'U', 0, nn_nevp )   ! First interpolation of coarse values 
     159      CALL agrif_interp_lim3( 'V', 0, nn_nevp ) 
    175160#endif 
    176161      ! 
    177162      !------------------------------------------------------------------------------! 
    178       ! 1) Ice strength (zpresh)                                ! 
    179       !------------------------------------------------------------------------------! 
     163      ! 0) mask at F points for the ice 
     164      !------------------------------------------------------------------------------! 
     165      ! ocean/land mask 
     166      DO jj = 1, jpjm1 
     167         DO ji = 1, jpim1      ! NO vector opt. 
     168            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 
     169         END DO 
     170      END DO 
     171      CALL lbc_lnk( zfmask, 'F', 1._wp ) 
     172 
     173      ! Lateral boundary conditions on velocity (modify zfmask) 
     174      zwf(:,:) = zfmask(:,:) 
     175      DO jj = 2, jpjm1 
     176         DO ji = fs_2, fs_jpim1   ! vector opt. 
     177            IF( zfmask(ji,jj) == 0._wp ) THEN 
     178               zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 
     179            ENDIF 
     180         END DO 
     181      END DO 
     182      DO jj = 2, jpjm1 
     183         IF( zfmask(1,jj) == 0._wp ) THEN 
     184            zfmask(1  ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 
     185         ENDIF 
     186         IF( zfmask(jpi,jj) == 0._wp ) THEN 
     187            zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
     188         ENDIF 
     189      END DO 
     190      DO ji = 2, jpim1 
     191         IF( zfmask(ji,1) == 0._wp ) THEN 
     192            zfmask(ji,1  ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 
     193         ENDIF 
     194         IF( zfmask(ji,jpj) == 0._wp ) THEN 
     195            zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 
     196         ENDIF 
     197      END DO 
     198      CALL lbc_lnk( zfmask, 'F', 1._wp ) 
     199 
     200      !------------------------------------------------------------------------------! 
     201      ! 1) define some variables and initialize arrays 
     202      !------------------------------------------------------------------------------! 
     203      zrhoco = rau0 * rn_cio  
     204 
     205      ! ecc2: square of yield ellipse eccenticrity 
     206      ecc2    = rn_ecc * rn_ecc 
     207      z1_ecc2 = 1._wp / ecc2 
     208 
     209      ! Time step for subcycling 
     210      zdtevp   = rdt_ice / REAL( nn_nevp ) 
     211      z1_dtevp = 1._wp / zdtevp 
     212 
     213      ! alpha parameters (Bouillon 2009) 
     214      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
     215      zalph2 = zalph1 * z1_ecc2 
     216 
     217      ! alpha and beta parameters (Bouillon 2013) 
     218      !!zalph1 = 40. 
     219      !!zalph2 = 40. 
     220      !!zbeta  = 3000. 
     221      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001) 
     222 
     223      z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     224      z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     225 
     226      ! Initialise stress tensor  
     227      zs1 (:,:) = stress1_i (:,:)  
     228      zs2 (:,:) = stress2_i (:,:) 
     229      zs12(:,:) = stress12_i(:,:) 
     230 
     231      ! Ice strength 
     232      CALL lim_itd_me_icestrength( nn_icestr ) 
     233 
     234      ! scale factors 
     235      DO jj = 2, jpjm1 
     236         DO ji = fs_2, fs_jpim1 
     237            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) ) 
     238            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) ) 
     239         END DO 
     240      END DO 
     241             
    180242      ! 
    181       ! Put every vector to 0 
    182       delta_i(:,:) = 0._wp   ; 
    183       zpresh (:,:) = 0._wp   ;   
    184       zpreshc(:,:) = 0._wp 
    185       u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp 
    186       divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp 
    187       shear_i(:,:) = 0._wp 
    188  
    189 #if defined key_lim3 
    190       CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points 
    191 #endif 
    192  
    193       DO jj = k_j1 , k_jpj       ! Ice mass and temp variables 
    194          DO ji = 1 , jpi 
    195 #if defined key_lim3 
    196             zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj) 
    197 #endif 
    198 #if defined key_lim2 
    199             zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 
    200 #endif 
    201             ! zmask = 1 where there is ice or on land 
    202             zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 
    203          END DO 
    204       END DO 
    205  
    206       ! Ice strength on grid cell corners (zpreshc) 
    207       ! needed for calculation of shear stress  
    208       DO jj = k_j1+1, k_jpj-1 
    209          DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 
    210             zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   & 
    211                &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1) 
    212             zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   & 
    213                &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     & 
    214                &             ) / MAX( zstms, zepsi ) 
    215          END DO 
    216       END DO 
    217       CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 
    218       ! 
    219243      !------------------------------------------------------------------------------! 
    220244      ! 2) Wind / ocean stress, mass terms, coriolis terms 
    221245      !------------------------------------------------------------------------------! 
    222       ! 
    223       !  Wind stress, coriolis and mass terms on the sides of the squares         
    224       !  zfrld1: lead fraction on U-points                                       
    225       !  zfrld2: lead fraction on V-points                                      
    226       !  zmass1: ice/snow mass on U-points                                     
    227       !  zmass2: ice/snow mass on V-points                                    
    228       !  zcorl1: Coriolis parameter on U-points                              
    229       !  zcorl2: Coriolis parameter on V-points                             
    230       !  (ztagnx,ztagny): wind stress on U/V points                        
    231       !  v_oce1: ocean v component on u points                           
    232       !  u_oce2: ocean u component on v points                          
    233246 
    234247      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
    … …  
    242255         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    243256         ! 
    244          zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:) ) * r1_rau0 
     257         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 
    245258         ! 
    246259      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
    … …  
    248261      ENDIF 
    249262 
    250       DO jj = k_j1+1, k_jpj-1 
     263      DO jj = 2, jpjm1 
    251264         DO ji = fs_2, fs_jpim1 
    252265 
    253             zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
    254             zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
    255             zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
    256  
    257             zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj) 
    258             zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 
    259             zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  ) 
    260             zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 
    261  
    262             ! Leads area. 
    263             zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi ) 
    264             zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi ) 
    265  
    266             ! Mass, coriolis coeff. and currents 
    267             zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 
    268             zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 
    269             zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   & 
    270                &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 
    271             zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   & 
    272                &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 
    273             ! 
    274             ! Ocean has no slip boundary condition 
    275             v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      & 
    276                &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  & 
    277                &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)   
    278  
    279             u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      & 
    280                &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  & 
    281                &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    282  
    283             ! Wind stress at U,V-point 
    284             ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj) 
    285             ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj) 
    286  
    287             ! Computation of the velocity field taking into account the ice internal interaction. 
    288             ! Terms that are independent of the velocity field. 
    289  
    290             ! SB On utilise maintenant le gradient de la pente de l'ocean 
    291             ! include it later 
    292  
    293             zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
    294             zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
    295  
    296             za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 
    297             za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy 
     266            ! ice fraction at U-V points 
     267            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) 
     268            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) 
     269 
     270            ! Ice/snow mass at U-V points 
     271            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
     272            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
     273            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
     274            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
     275            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     276 
     277            ! Ocean currents at U-V points 
     278            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    & 
     279               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     280             
     281            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    & 
     282               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     283 
     284            ! Coriolis at T points (m*f) 
     285            zmf(ji,jj)      = zm1 * ff_t(ji,jj) 
     286 
     287            ! m/dt 
     288            zmU_t(ji,jj)    = zmassU * z1_dtevp 
     289            zmV_t(ji,jj)    = zmassV * z1_dtevp 
     290 
     291            ! Drag ice-atm. 
     292            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     293            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
     294 
     295            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
     296            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
     297            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
     298 
     299            ! masks 
     300            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
     301            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
     302 
     303            ! switches 
     304            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 
     305            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 
    298306 
    299307         END DO 
    300308      END DO 
    301  
     309      CALL lbc_lnk( zmf, 'T', 1. ) 
    302310      ! 
    303311      !------------------------------------------------------------------------------! 
    … …  
    305313      !------------------------------------------------------------------------------! 
    306314      ! 
    307       ! Time step for subcycling 
    308       dtevp  = rdt_ice / nn_nevp 
    309 #if defined key_lim3 
    310       dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 
    311 #else 
    312       dtotel = dtevp / ( 2._wp * telast ) 
    313 #endif 
    314       z1_dtotel = 1._wp / ( 1._wp + dtotel ) 
    315       z1_dtevp  = 1._wp / dtevp 
    316       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 
    317       ecc2 = rn_ecc * rn_ecc 
    318       ecci = 1. / ecc2 
    319  
    320       !-Initialise stress tensor  
    321       zs1 (:,:) = stress1_i (:,:)  
    322       zs2 (:,:) = stress2_i (:,:) 
    323       zs12(:,:) = stress12_i(:,:) 
    324  
    325315      !                                               !----------------------! 
    326316      DO jter = 1 , nn_nevp                           !    loop over jter    ! 
    327317         !                                            !----------------------!         
    328          DO jj = k_j1, k_jpj-1 
    329             zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step 
    330             zv_ice(:,jj) = v_ice(:,jj) 
    331          END DO 
    332  
    333          DO jj = k_j1+1, k_jpj-1 
    334             DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask 
    335  
    336                !   
    337                !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 
    338                !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 
    339                !- zds(:,:): shear on northeast corner of grid cells 
    340                ! 
    341                !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,  
    342                !                      there are many repeated calculations.  
    343                !                      Speed could be improved by regrouping terms. For 
    344                !                      the moment, however, the stress is on clarity of coding to avoid 
    345                !                      bugs (Martin, for Miguel). 
    346                ! 
    347                !- ALSO: arrays zdt, zds and delta could  
    348                !  be removed in the future to minimise memory demand. 
    349                ! 
    350                !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of 
    351                !              grid cells, exactly as in the B grid case. For simplicity, the indexation on 
    352                !              the corners is the same as in the B grid. 
    353                ! 
    354                ! 
    355                divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    356                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    357                   &            ) * r1_e1e2t(ji,jj) 
    358  
    359                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    360                   &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    361                   &         ) * r1_e1e2t(ji,jj) 
    362  
    363                ! 
     318         IF(ln_ctl) THEN   ! Convergence test 
     319            DO jj = 1, jpjm1 
     320               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
     321               zv_ice(:,jj) = v_ice(:,jj) 
     322            END DO 
     323         ENDIF 
     324 
     325         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
     326         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 
     327            DO ji = 1, jpim1 
     328 
     329               ! shear at F points 
    364330               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)   & 
    365331                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    366                   &         ) * r1_e1e2f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   & 
    367                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
    368  
    369  
    370                v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    371                   &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   & 
    372                   &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)  
    373  
    374                u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    375                   &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   & 
    376                   &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    377             END DO 
    378          END DO 
    379  
    380          CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond. 
     332                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
     333 
     334            END DO 
     335         END DO 
     336         CALL lbc_lnk( zds, 'F', 1. ) 
     337 
     338         DO jj = 2, jpjm1 
     339            DO ji = 2, jpim1 ! no vector loop 
     340 
     341               ! shear**2 at T points (doc eq. A16) 
     342               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     343                  &   + 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)  & 
     344                  &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
     345               
     346               ! divergence at T points 
     347               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     348                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     349                  &    ) * r1_e1e2t(ji,jj) 
     350               zdiv2 = zdiv * zdiv 
     351                
     352               ! tension at T points 
     353               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)   & 
     354                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     355                  &   ) * r1_e1e2t(ji,jj) 
     356               zdt2 = zdt * zdt 
     357                
     358               ! delta at T points 
     359               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     360 
     361               ! P/delta at T points 
     362               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
     363                
     364               ! stress at T points 
     365               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 
     366               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 
     367              
     368            END DO 
     369         END DO 
     370         CALL lbc_lnk( zp_delt, 'T', 1. ) 
     371 
     372         DO jj = 1, jpjm1 
     373            DO ji = 1, jpim1 
     374 
     375               ! P/delta at F points 
     376               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) ) 
     377                
     378               ! stress at F points 
     379               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 
     380 
     381            END DO 
     382         END DO 
     383         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
     384  
     385 
     386         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
     387         DO jj = 2, jpjm1 
     388            DO ji = fs_2, fs_jpim1                
     389 
     390               ! U points 
     391               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
     392                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     393                  &                    ) * r1_e2u(ji,jj)                                                                      & 
     394                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
     395                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
     396                  &                  ) * r1_e1e2u(ji,jj) 
     397 
     398               ! V points 
     399               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
     400                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     401                  &                    ) * r1_e1v(ji,jj)                                                                      & 
     402                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
     403                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
     404                  &                  ) * r1_e1e2v(ji,jj) 
     405 
     406               ! u_ice at V point 
     407               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     408                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     409                
     410               ! v_ice at U point 
     411               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
     412                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     413 
     414            END DO 
     415         END DO 
     416         ! 
     417         ! --- Computation of ice velocity --- ! 
     418         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 
     419         !  Bouillon et al. 2009 (eq 34-35) => stable 
     420         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations 
     421             
     422            DO jj = 2, jpjm1 
     423               DO ji = fs_2, fs_jpim1 
     424 
     425                  ! tau_io/(v_oce - v_ice) 
     426                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     427                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     428 
     429                  ! tau_bottom/v_ice 
     430                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     431                  zTauB = - tau_icebfr(ji,jj) / zvel 
     432 
     433                  ! Coriolis at V-points (energy conserving formulation) 
     434                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     435                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     436                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     437 
     438                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     439                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     440 
     441                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     442                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     443                   
     444                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     445                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
     446                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     447                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     448                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     449                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
     450                     &           ) * zmaskV(ji,jj) 
     451                  ! Bouillon 2013 
     452                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     453                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
     454                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
     455 
     456               END DO 
     457            END DO 
     458            CALL lbc_lnk( v_ice, 'V', -1. ) 
     459             
     460#if defined key_agrif 
     461!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp ) 
     462            CALL agrif_interp_lim3( 'V' ) 
     463#endif 
     464            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' ) 
     465 
     466            DO jj = 2, jpjm1 
     467               DO ji = fs_2, fs_jpim1 
     468                                
     469                  ! tau_io/(u_oce - u_ice) 
     470                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     471                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     472 
     473                  ! tau_bottom/u_ice 
     474                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     475                  zTauB = - tau_icebfr(ji,jj) / zvel 
     476 
     477                  ! Coriolis at U-points (energy conserving formulation) 
     478                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     479                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     480                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     481                   
     482                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     483                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     484 
     485                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     486                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     487 
     488                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     489                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
     490                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     491                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     492                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     493                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin  
     494                     &           ) * zmaskU(ji,jj) 
     495                  ! Bouillon 2013 
     496                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     497                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
     498                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
     499               END DO 
     500            END DO 
     501            CALL lbc_lnk( u_ice, 'U', -1. ) 
     502             
     503#if defined key_agrif 
     504!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp ) 
     505            CALL agrif_interp_lim3( 'U' ) 
     506#endif 
     507            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' ) 
     508 
     509         ELSE ! odd iterations 
     510 
     511            DO jj = 2, jpjm1 
     512               DO ji = fs_2, fs_jpim1 
     513 
     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 
     518                  ! tau_bottom/u_ice 
     519                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
     520                  zTauB = - tau_icebfr(ji,jj) / zvel 
     521 
     522                  ! Coriolis at U-points (energy conserving formulation) 
     523                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     524                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     525                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     526 
     527                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     528                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     529 
     530                  ! landfast switch => 0 = static friction ; 1 = sliding friction 
     531                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     532 
     533                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     534                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity 
     535                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     536                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     537                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     538                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
     539                     &           ) * zmaskU(ji,jj) 
     540                  ! Bouillon 2013 
     541                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
     542                  !!   &           + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
     543                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
     544               END DO 
     545            END DO 
     546            CALL lbc_lnk( u_ice, 'U', -1. ) 
     547             
     548#if defined key_agrif 
     549!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp ) 
     550            CALL agrif_interp_lim3( 'U' ) 
     551#endif 
     552            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' ) 
     553 
     554            DO jj = 2, jpjm1 
     555               DO ji = fs_2, fs_jpim1 
    381556          
    382          DO jj = k_j1+1, k_jpj-1 
    383             DO ji = fs_2, fs_jpim1 
    384  
    385                !- Calculate Delta at centre of grid cells 
    386                zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   & 
    387                   &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   & 
    388                   &            ) * r1_e1e2t(ji,jj) 
    389  
    390                delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    391                delta_i(ji,jj) = delta + rn_creepl 
    392  
    393                !- Calculate Delta on corners 
    394                zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    395                   &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    396                   &    ) * r1_e1e2f(ji,jj) 
    397  
    398                zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    399                   &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    400                   &    ) * r1_e1e2f(ji,jj) 
    401  
    402                zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 
    403  
    404                !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 
    405                zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
    406                   &          ) * z1_dtotel 
    407                zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
    408                   &          ) * z1_dtotel 
    409                !-Calculate stress tensor component zs12 at corners 
    410                zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  & 
    411                   &          ) * z1_dtotel  
    412  
    413             END DO 
    414          END DO 
    415  
    416          CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    417   
    418          ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
    419          DO jj = k_j1+1, k_jpj-1 
    420             DO ji = fs_2, fs_jpim1 
    421                !- contribution of zs1, zs2 and zs12 to zf1 
    422                zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  & 
    423                   &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          & 
    424                   &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  & 
    425                   &                ) * r1_e1e2u(ji,jj) 
    426                ! contribution of zs1, zs2 and zs12 to zf2 
    427                zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  & 
    428                   &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          & 
    429                   &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  & 
    430                   &               )  * r1_e1e2v(ji,jj) 
    431             END DO 
    432          END DO 
    433          ! 
    434          ! Computation of ice velocity 
    435          ! 
    436          ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 
    437          ! 
    438          IF (MOD(jter,2).eq.0) THEN  
    439  
    440             DO jj = k_j1+1, k_jpj-1 
    441                DO ji = fs_2, fs_jpim1 
    442                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 
    443                   z0           = zmass1(ji,jj) * z1_dtevp 
    444  
    445                   ! SB modif because ocean has no slip boundary condition 
    446                   zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     & 
    447                      &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
    448                      &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    449                   za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
    450                      &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
    451                   zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
    452                   zcca         = z0 + za 
    453                   zccb         = zcorl1(ji,jj) 
    454                   u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
     557                  ! tau_io/(v_oce - v_ice) 
     558                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     559                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     560 
     561                  ! tau_bottom/v_ice 
     562                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
     563                  ztauB = - tau_icebfr(ji,jj) / zvel 
     564                   
     565                  ! Coriolis at V-points (energy conserving formulation) 
     566                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     567                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     568                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     569 
     570                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     571                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     572 
     573                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction 
     574                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     575                   
     576                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     577                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity 
     578                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     579                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast 
     580                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
     581                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
     582                     &           ) * zmaskV(ji,jj) 
     583                  ! Bouillon 2013 
     584                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
     585                  !!   &           + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
     586                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    455587               END DO 
    456588            END DO 
    457  
    458             CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    459 #if defined key_agrif && defined key_lim2 
    460             CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
     589            CALL lbc_lnk( v_ice, 'V', -1. ) 
     590             
     591#if defined key_agrif 
     592!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp ) 
     593            CALL agrif_interp_lim3( 'V' ) 
    461594#endif 
    462 #if defined key_bdy 
    463          CALL bdy_ice_lim_dyn( 'U' ) 
    464 #endif          
    465  
    466             DO jj = k_j1+1, k_jpj-1 
    467                DO ji = fs_2, fs_jpim1 
    468  
    469                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
    470                   z0           = zmass2(ji,jj) * z1_dtevp 
    471                   ! SB modif because ocean has no slip boundary condition 
    472                   zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
    473                      &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
    474                      &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    475                   za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &  
    476                      &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
    477                   zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
    478                   zcca         = z0 + za 
    479                   zccb         = zcorl2(ji,jj) 
    480                   v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    481                END DO 
    482             END DO 
    483  
    484             CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    485 #if defined key_agrif && defined key_lim2 
    486             CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    487 #endif 
    488 #if defined key_bdy 
    489          CALL bdy_ice_lim_dyn( 'V' ) 
    490 #endif          
    491  
    492          ELSE  
    493             DO jj = k_j1+1, k_jpj-1 
    494                DO ji = fs_2, fs_jpim1 
    495                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
    496                   z0           = zmass2(ji,jj) * z1_dtevp 
    497                   ! SB modif because ocean has no slip boundary condition 
    498                   zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
    499                      &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
    500                      &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)    
    501  
    502                   za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
    503                      &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
    504                   zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
    505                   zcca         = z0 + za 
    506                   zccb         = zcorl2(ji,jj) 
    507                   v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    508                END DO 
    509             END DO 
    510  
    511             CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    512 #if defined key_agrif && defined key_lim2 
    513             CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    514 #endif 
    515 #if defined key_bdy 
    516          CALL bdy_ice_lim_dyn( 'V' ) 
    517 #endif          
    518  
    519             DO jj = k_j1+1, k_jpj-1 
    520                DO ji = fs_2, fs_jpim1 
    521                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 
    522                   z0           = zmass1(ji,jj) * z1_dtevp 
    523                   zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       & 
    524                      &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
    525                      &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    526  
    527                   za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
    528                      &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
    529                   zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
    530                   zcca         = z0 + za 
    531                   zccb         = zcorl1(ji,jj) 
    532                   u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
    533                END DO 
    534             END DO 
    535  
    536             CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    537 #if defined key_agrif && defined key_lim2 
    538             CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
    539 #endif 
    540 #if defined key_bdy 
    541          CALL bdy_ice_lim_dyn( 'U' ) 
    542 #endif          
     595            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' ) 
    543596 
    544597         ENDIF 
    545598          
    546          IF(ln_ctl) THEN 
    547             !---  Convergence test. 
    548             DO jj = k_j1+1 , k_jpj-1 
     599         IF(ln_ctl) THEN   ! Convergence test 
     600            DO jj = 2 , jpjm1 
    549601               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    550602            END DO 
    551             zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) ) 
     603            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 
    552604            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
    553605         ENDIF 
    554  
     606         ! 
    555607         !                                                ! ==================== ! 
    556608      END DO                                              !  end loop over jter  ! 
    … …  
    558610      ! 
    559611      !------------------------------------------------------------------------------! 
    560       ! 4) Prevent ice velocities when the ice is thin 
    561       !------------------------------------------------------------------------------! 
    562       ! If the ice volume is below zvmin then ice velocity should equal the 
    563       ! ocean velocity. This prevents high velocity when ice is thin 
    564       DO jj = k_j1+1, k_jpj-1 
    565          DO ji = fs_2, fs_jpim1 
    566             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    567                u_ice(ji,jj) = u_oce(ji,jj) 
    568                v_ice(ji,jj) = v_oce(ji,jj) 
    569             ENDIF 
     612      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)  
     613      !------------------------------------------------------------------------------! 
     614      DO jj = 1, jpjm1 
     615         DO ji = 1, jpim1 
     616 
     617            ! shear at F points 
     618            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)   & 
     619               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     620               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
     621 
     622         END DO 
     623      END DO            
     624      CALL lbc_lnk( zds, 'F', 1. ) 
     625       
     626      DO jj = 2, jpjm1 
     627         DO ji = 2, jpim1 ! no vector loop 
     628             
     629            ! tension**2 at T points 
     630            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)   & 
     631               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     632               &   ) * r1_e1e2t(ji,jj) 
     633            zdt2 = zdt * zdt 
     634             
     635            ! shear**2 at T points (doc eq. A16) 
     636            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     637               &   + 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)  & 
     638               &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
     639             
     640            ! shear at T points 
     641            shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     642 
     643            ! divergence at T points 
     644            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     645               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     646               &            ) * r1_e1e2t(ji,jj) 
     647             
     648            ! delta at T points 
     649            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
     650            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
     651            delta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     652 
    570653         END DO 
    571654      END DO 
    572  
    573       CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 
    574  
    575 #if defined key_agrif && defined key_lim2 
    576       CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 
    577       CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 
    578 #endif 
    579 #if defined key_bdy 
    580       CALL bdy_ice_lim_dyn( 'U' ) 
    581       CALL bdy_ice_lim_dyn( 'V' ) 
    582 #endif          
    583  
    584       DO jj = k_j1+1, k_jpj-1  
    585          DO ji = fs_2, fs_jpim1 
    586             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    587                v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     & 
    588                   &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   & 
    589                   &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    590  
    591                u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    592                   &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   & 
    593                   &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    594             ENDIF  
    595          END DO 
    596       END DO 
    597  
    598       CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 
    599  
    600       ! Recompute delta, shear and div, inputs for mechanical redistribution  
    601       DO jj = k_j1+1, k_jpj-1 
    602          DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask 
    603             !- divu_i(:,:), zdt(:,:): divergence and tension at centre  
    604             !- zds(:,:): shear on northeast corner of grid cells 
    605             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    606  
    607                divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   & 
    608                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   & 
    609                   &            ) * r1_e1e2t(ji,jj) 
    610  
    611                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  & 
    612                   &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  & 
    613                   &         ) * r1_e1e2t(ji,jj) 
    614                ! 
    615                ! SB modif because ocean has no slip boundary condition  
    616                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)  & 
    617                   &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    618                   &         ) * r1_e1e2f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     & 
    619                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
    620  
    621                zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    & 
    622                   &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e1e2t(ji,jj) 
    623  
    624                delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    625                delta_i(ji,jj) = delta + rn_creepl 
    626              
    627             ENDIF 
    628          END DO 
    629       END DO 
    630       ! 
    631       !------------------------------------------------------------------------------! 
    632       ! 5) Store stress tensor and its invariants 
    633       !------------------------------------------------------------------------------! 
    634       ! * Invariants of the stress tensor are required for limitd_me 
    635       !   (accelerates convergence and improves stability) 
    636       DO jj = k_j1+1, k_jpj-1 
    637          DO ji = fs_2, fs_jpim1 
    638             zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &    
    639                &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e1e2t(ji,jj)  
    640             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 
    641          END DO 
    642       END DO 
    643  
    644       ! Lateral boundary condition 
    645       CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. ) 
    646  
    647       ! * Store the stress tensor for the next time step 
     655      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 
     656       
     657      ! --- Store the stress tensor for the next time step --- ! 
    648658      stress1_i (:,:) = zs1 (:,:) 
    649659      stress2_i (:,:) = zs2 (:,:) 
    650660      stress12_i(:,:) = zs12(:,:) 
    651  
    652661      ! 
    653       !------------------------------------------------------------------------------! 
    654       ! 6) Control prints of residual and charge ellipse 
     662 
     663      !------------------------------------------------------------------------------! 
     664      ! 5) Control prints of residual and charge ellipse 
    655665      !------------------------------------------------------------------------------! 
    656666      ! 
    … …  
    672682            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. ' 
    673683            CALL prt_ctl_info(charout) 
    674             DO jj = k_j1+1, k_jpj-1 
     684            DO jj = 2, jpjm1 
    675685               DO ji = 2, jpim1 
    676                   IF (zpresh(ji,jj) > 1.0) THEN 
    677                      sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )  
    678                      sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
     686                  IF (strength(ji,jj) > 1.0) THEN 
     687                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )  
     688                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
    679689                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 
    680690                     CALL prt_ctl_info(charout) 
    … …  
    686696         ENDIF 
    687697      ENDIF 
    688       ! 
    689       CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    690       CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    691       CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    692       CALL wrk_dealloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 ) 
     698      !      
     699       
     700      CALL wrk_dealloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt ) 
     701      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 
     702      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
     703      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     704      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
    693705 
    694706   END SUBROUTINE lim_rhg 
    … …  
    699711   !!---------------------------------------------------------------------- 
    700712CONTAINS 
    701    SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine 
    702       WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2 
     713   SUBROUTINE lim_rhg         ! Dummy routine 
     714      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?' 
    703715   END SUBROUTINE lim_rhg 
    704716#endif 
Note: See TracChangeset for help on using the changeset viewer.