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 5600 for branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 – NEMO

Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (9 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r5034 r5600  
    102102      !!                 and charge ellipse. 
    103103      !!                 The user should make sure that the parameters 
    104       !!                 nevp, telast and creepl maintain stress state 
     104      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state 
    105105      !!                 on the charge ellipse for plastic flow 
    106106      !!                 e.g. in the Canadian Archipelago 
     
    108108      !! References : Hunke and Dukowicz, JPO97 
    109109      !!              Bouillon et al., Ocean Modelling 2009 
    110       !!              Vancoppenolle et al., Ocean Modelling 2008 
    111110      !!------------------------------------------------------------------- 
    112111      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation 
     
    117116      CHARACTER (len=50) ::   charout 
    118117      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         ! 
    119       REAL(wp) ::   za, zstms, zmask   ! local scalars 
    120       REAL(wp) ::   zc1, zc2, zc3             ! ice mass 
    121  
    122       REAL(wp) ::   dtevp              ! time step for subcycling 
    123       REAL(wp) ::   dtotel, ecc2, ecci ! square of yield ellipse eccenticity 
    124       REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars 
    125       REAL(wp) ::   zu_ice2, zv_ice1   ! 
    126       REAL(wp) ::   zddc, zdtc         ! delta on corners and on centre 
    127       REAL(wp) ::   zdst               ! shear at the center of the grid point 
    128       REAL(wp) ::   zdsshx, zdsshy     ! term for the gradient of ocean surface 
    129       REAL(wp) ::   sigma1, sigma2     ! internal ice stress 
     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 
    130129 
    131130      REAL(wp) ::   zresm         ! Maximal error on ice velocity 
    132       REAL(wp) ::   zdummy        ! dummy argument 
    133131      REAL(wp) ::   zintb, zintn  ! dummy argument 
    134132 
     
    139137      REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points 
    140138      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays 
    141       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce1, v_oce1   ! ocean u/v component on U points                            
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2, v_oce2   ! ocean u/v component on V points 
     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 
    143141      REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point 
    144142      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses 
     143      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points 
    145144       
    146145      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells 
     
    152151                                                              !   ocean surface (ssh_m) if ice is not embedded 
    153152                                                              !   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 
    154156      !!------------------------------------------------------------------- 
    155157 
    156158      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    157       CALL wrk_alloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                ) 
     159      CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    158160      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    159161      CALL wrk_alloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
     
    161163#if  defined key_lim2 && ! defined key_lim2_vp 
    162164# if defined key_agrif 
    163      USE ice_2, vt_s => hsnm 
    164      USE ice_2, vt_i => hicm 
     165      USE ice_2, vt_s => hsnm 
     166      USE ice_2, vt_i => hicm 
    165167# else 
    166      vt_s => hsnm 
    167      vt_i => hicm 
     168      vt_s => hsnm 
     169      vt_i => hicm 
    168170# endif 
    169      at_i(:,:) = 1. - frld(:,:) 
     171      at_i(:,:) = 1. - frld(:,:) 
    170172#endif 
    171173#if defined key_agrif && defined key_lim2  
    172     CALL agrif_rhg_lim2_load      ! First interpolation of coarse values 
     174      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values 
    173175#endif 
    174176      ! 
     
    186188 
    187189#if defined key_lim3 
    188       CALL lim_itd_me_icestrength( ridge_scheme_swi )      ! LIM-3: Ice strength on T-points 
    189 #endif 
    190  
    191 !CDIR NOVERRCHK 
     190      CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points 
     191#endif 
     192 
    192193      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables 
    193 !CDIR NOVERRCHK 
    194194         DO ji = 1 , jpi 
    195195#if defined key_lim3 
    196             zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) 
     196            zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj) 
    197197#endif 
    198198#if defined key_lim2 
    199             zpresh(ji,jj) = tms(ji,jj) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 
    200 #endif 
    201             ! tmi = 1 where there is ice or on land 
    202             tmi(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * tms(ji,jj) 
     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) 
    203203         END DO 
    204204      END DO 
     
    206206      ! Ice strength on grid cell corners (zpreshc) 
    207207      ! needed for calculation of shear stress  
    208 !CDIR NOVERRCHK 
    209208      DO jj = k_j1+1, k_jpj-1 
    210 !CDIR NOVERRCHK 
    211209         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 
    212             zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 
    213                &              tms(ji,jj+1)   * wght(ji+1,jj+1,1,2) + & 
    214                &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
    215                &              tms(ji,jj)     * wght(ji+1,jj+1,1,1) 
    216             zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 
    217                &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + & 
    218                &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &  
    219                &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   & 
    220                &             ) / MAX( zstms, epsd ) 
     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 ) 
    221215         END DO 
    222216      END DO 
    223  
    224217      CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 
    225218      ! 
     
    236229      !  zcorl2: Coriolis parameter on V-points                             
    237230      !  (ztagnx,ztagny): wind stress on U/V points                        
    238       !  u_oce1: ocean u component on u points                            
    239231      !  v_oce1: ocean v component on u points                           
    240232      !  u_oce2: ocean u component on v points                          
    241       !  v_oce2: ocean v component on v points                         
    242233 
    243234      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
    244           !                                             
    245           ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
    246           !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 
     235         !                                             
     236         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
     237         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 
    247238         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp      
    248           ! 
    249           ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
    250           !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
     239         ! 
     240         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
     241         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
    251242         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    252           ! 
     243         ! 
    253244         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0 
    254           ! 
     245         ! 
    255246      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
    256247         zpice(:,:) = ssh_m(:,:) 
     
    260251         DO ji = fs_2, fs_jpim1 
    261252 
    262             zc1 = tms(ji  ,jj  ) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
    263             zc2 = tms(ji+1,jj  ) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
    264             zc3 = tms(ji  ,jj+1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
    265  
    266             zt11 = tms(ji  ,jj) * e1t(ji  ,jj) 
    267             zt12 = tms(ji+1,jj) * e1t(ji+1,jj) 
    268             zt21 = tms(ji,jj  ) * e2t(ji,jj  ) 
    269             zt22 = tms(ji,jj+1) * e2t(ji,jj+1) 
     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) 
    270261 
    271262            ! Leads area. 
    272             zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd ) 
    273             zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd ) 
     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 ) 
    274265 
    275266            ! Mass, coriolis coeff. and currents 
    276             zmass1(ji,jj) = ( zt12*zc1 + zt11*zc2 ) / (zt11+zt12+epsd) 
    277             zmass2(ji,jj) = ( zt22*zc1 + zt21*zc3 ) / (zt21+zt22+epsd) 
    278             zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) )   & 
    279                &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd ) 
    280             zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) )   & 
    281                &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd ) 
     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 ) 
    282273            ! 
    283             u_oce1(ji,jj)  = u_oce(ji,jj) 
    284             v_oce2(ji,jj)  = v_oce(ji,jj) 
    285  
    286274            ! Ocean has no slip boundary condition 
    287             v_oce1(ji,jj)  = 0.5*( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)    & 
    288                &                 +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) & 
    289                &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj 
    290  
    291             u_oce2(ji,jj)  = 0.5*((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)     & 
    292                &                 +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) & 
    293                &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
     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) 
    294282 
    295283            ! Wind stress at U,V-point 
     
    303291            ! include it later 
    304292 
    305             zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
    306             zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
     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) 
    307295 
    308296            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 
     
    318306      ! 
    319307      ! Time step for subcycling 
    320       dtevp  = rdt_ice / nevp 
     308      dtevp  = rdt_ice / nn_nevp 
     309#if defined key_lim3 
     310      dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 
     311#else 
    321312      dtotel = dtevp / ( 2._wp * telast ) 
    322  
     313#endif 
     314      z1_dtotel = 1._wp / ( 1._wp + dtotel ) 
     315      z1_dtevp  = 1._wp / dtevp 
    323316      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 
    324       ecc2 = ecc * ecc 
     317      ecc2 = rn_ecc * rn_ecc 
    325318      ecci = 1. / ecc2 
    326319 
     
    331324 
    332325      !                                               !----------------------! 
    333       DO jter = 1 , nevp                              !    loop over jter    ! 
     326      DO jter = 1 , nn_nevp                           !    loop over jter    ! 
    334327         !                                            !----------------------!         
    335328         DO jj = k_j1, k_jpj-1 
     
    339332 
    340333         DO jj = k_j1+1, k_jpj-1 
    341             DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi 
     334            DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask 
    342335 
    343336               !   
     
    360353               ! 
    361354               ! 
    362                divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
    363                   &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  & 
    364                   &             +e1v(ji,jj)*v_ice(ji,jj)                      & 
    365                   &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  & 
    366                   &             )                                             & 
    367                   &            / area(ji,jj) 
    368  
    369                zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    & 
    370                   &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                & 
    371                   &           )*e2t(ji,jj)*e2t(ji,jj)                      & 
    372                   &          -( v_ice(ji,jj)/e1v(ji,jj)                    & 
    373                   &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                & 
    374                   &           )*e1t(ji,jj)*e1t(ji,jj)                      & 
    375                   &         )                                              & 
    376                   &        / area(ji,jj) 
     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_e12t(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_e12t(ji,jj) 
    377362 
    378363               ! 
    379                zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1)                & 
    380                   &            -u_ice(ji,jj)/e1u(ji,jj)                    & 
    381                   &           )*e1f(ji,jj)*e1f(ji,jj)                      & 
    382                   &          +( v_ice(ji+1,jj)/e2v(ji+1,jj)                & 
    383                   &            -v_ice(ji,jj)/e2v(ji,jj)                    & 
    384                   &           )*e2f(ji,jj)*e2f(ji,jj)                      & 
    385                   &         )                                              & 
    386                   &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 
    387                   &        * tmi(ji,jj) * tmi(ji,jj+1)                     & 
    388                   &        * tmi(ji+1,jj) * tmi(ji+1,jj+1) 
    389  
    390  
    391                v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   & 
    392                   &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 
    393                   &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)  
    394  
    395                u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   & 
    396                   &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 
    397                   &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
    398  
    399             END DO 
    400          END DO 
    401          CALL lbc_lnk( v_ice1, 'U', -1. )   ;   CALL lbc_lnk( u_ice2, 'V', -1. )      ! lateral boundary cond. 
    402  
    403 !CDIR NOVERRCHK 
     364               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)   & 
     365                  &         + ( 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_e12f(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. 
     381          
    404382         DO jj = k_j1+1, k_jpj-1 
    405 !CDIR NOVERRCHK 
    406383            DO ji = fs_2, fs_jpim1 
    407384 
    408385               !- Calculate Delta at centre of grid cells 
    409                zdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          & 
    410                   &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          & 
    411                   &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          & 
    412                   &          - e1v(ji, jj-1) * u_ice2(ji,jj-1)          & 
    413                   &          )                                          & 
    414                   &         / area(ji,jj) 
    415  
    416                delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )   
    417                delta_i(ji,jj) = delta + creepl 
    418                !-Calculate stress tensor components zs1 and zs2  
    419                !-at centre of grid cells (see section 3.5 of CICE user's guide). 
    420                zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) )  & 
    421                   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 
    422                zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) )  & 
    423                   &         / ( 1._wp + dtotel ) 
    424  
    425             END DO 
    426          END DO 
    427  
    428          CALL lbc_lnk( zs1(:,:), 'T', 1. ) 
    429          CALL lbc_lnk( zs2(:,:), 'T', 1. ) 
    430  
    431 !CDIR NOVERRCHK 
    432          DO jj = k_j1+1, k_jpj-1 
    433 !CDIR NOVERRCHK 
    434             DO ji = fs_2, fs_jpim1 
     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_e12t(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 
    435393               !- Calculate Delta on corners 
    436                zddc  =      ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1)                & 
    437                   &            -v_ice1(ji,jj)/e1u(ji,jj)                    & 
    438                   &           )*e1f(ji,jj)*e1f(ji,jj)                       & 
    439                   &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                & 
    440                   &            -u_ice2(ji,jj)/e2v(ji,jj)                    & 
    441                   &           )*e2f(ji,jj)*e2f(ji,jj)                       & 
    442                   &         )                                               & 
    443                   &        / ( e1f(ji,jj) * e2f(ji,jj) ) 
    444  
    445                zdtc  =      (-( v_ice1(ji,jj+1)/e1u(ji,jj+1)                & 
    446                   &            -v_ice1(ji,jj)/e1u(ji,jj)                    & 
    447                   &           )*e1f(ji,jj)*e1f(ji,jj)                       & 
    448                   &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                & 
    449                   &            -u_ice2(ji,jj)/e2v(ji,jj)                    & 
    450                   &           )*e2f(ji,jj)*e2f(ji,jj)                       & 
    451                   &         )                                               & 
    452                   &        / ( e1f(ji,jj) * e2f(ji,jj) ) 
    453  
    454                zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl 
    455  
    456                !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 
    457                zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  & 
    458                   &          ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) )  & 
    459                   &          / ( 1.0 + dtotel )  
    460  
    461             END DO ! ji 
    462          END DO ! jj 
    463  
    464          CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
    465  
     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_e12f(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_e12f(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  
    466418         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
    467419         DO jj = k_j1+1, k_jpj-1 
    468420            DO ji = fs_2, fs_jpim1 
    469421               !- contribution of zs1, zs2 and zs12 to zf1 
    470                zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) & 
    471                   &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) & 
    472                   &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) & 
    473                   &             ) / ( e1u(ji,jj)*e2u(ji,jj) ) 
     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_e12u(ji,jj) 
    474426               ! contribution of zs1, zs2 and zs12 to zf2 
    475                zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) & 
    476                   &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) & 
    477                   &              + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 -    & 
    478                   zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) & 
    479                   &             ) / ( e1v(ji,jj)*e2v(ji,jj) ) 
     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_e12v(ji,jj) 
    480431            END DO 
    481432         END DO 
     
    487438         IF (MOD(jter,2).eq.0) THEN  
    488439 
    489 !CDIR NOVERRCHK 
    490440            DO jj = k_j1+1, k_jpj-1 
    491 !CDIR NOVERRCHK 
    492441               DO ji = fs_2, fs_jpim1 
    493                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 
    494                   z0           = zmass1(ji,jj)/dtevp 
     442                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 
     443                  z0           = zmass1(ji,jj) * z1_dtevp 
    495444 
    496445                  ! SB modif because ocean has no slip boundary condition 
    497                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)         & 
    498                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   & 
    499                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
    500                   za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 
    501                      (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj)) 
    502                   zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 
    503                      za*(u_oce1(ji,jj)) 
    504                   zcca         = z0+za 
     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 
    505453                  zccb         = zcorl1(ji,jj) 
    506                   u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask  
    507  
     454                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
    508455               END DO 
    509456            END DO 
     
    511458            CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    512459#if defined key_agrif && defined key_lim2 
    513             CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
     460            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
    514461#endif 
    515462#if defined key_bdy 
     
    517464#endif          
    518465 
    519 !CDIR NOVERRCHK 
    520466            DO jj = k_j1+1, k_jpj-1 
    521 !CDIR NOVERRCHK 
    522467               DO ji = fs_2, fs_jpim1 
    523468 
    524                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 
    525                   z0           = zmass2(ji,jj)/dtevp 
     469                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
     470                  z0           = zmass2(ji,jj) * z1_dtevp 
    526471                  ! SB modif because ocean has no slip boundary condition 
    527                   zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)     & 
    528                      &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   & 
    529                      &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
    530                   za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + &  
    531                      (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 
    532                   zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + & 
    533                      za2ct(ji,jj) + za*(v_oce2(ji,jj)) 
    534                   zcca         = z0+za 
     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 
    535479                  zccb         = zcorl2(ji,jj) 
    536                   v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 
    537  
     480                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    538481               END DO 
    539482            END DO 
     
    541484            CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    542485#if defined key_agrif && defined key_lim2 
    543             CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
     486            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    544487#endif 
    545488#if defined key_bdy 
     
    548491 
    549492         ELSE  
    550 !CDIR NOVERRCHK 
    551493            DO jj = k_j1+1, k_jpj-1 
    552 !CDIR NOVERRCHK 
    553494               DO ji = fs_2, fs_jpim1 
    554                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 
    555                   z0           = zmass2(ji,jj)/dtevp 
     495                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
     496                  z0           = zmass2(ji,jj) * z1_dtevp 
    556497                  ! SB modif because ocean has no slip boundary condition 
    557                   zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)      & 
    558                      &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   & 
    559                      &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)    
    560  
    561                   za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 
    562                      (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj)) 
    563                   zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + & 
    564                      za2ct(ji,jj) + za*(v_oce2(ji,jj)) 
    565                   zcca         = z0+za 
     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 
    566506                  zccb         = zcorl2(ji,jj) 
    567                   v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+epsd)*zmask 
    568  
     507                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    569508               END DO 
    570509            END DO 
     
    572511            CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    573512#if defined key_agrif && defined key_lim2 
    574             CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
     513            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    575514#endif 
    576515#if defined key_bdy 
     
    578517#endif          
    579518 
    580 !CDIR NOVERRCHK 
    581519            DO jj = k_j1+1, k_jpj-1 
    582 !CDIR NOVERRCHK 
    583520               DO ji = fs_2, fs_jpim1 
    584                   zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 
    585                   z0           = zmass1(ji,jj)/dtevp 
    586                   zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      & 
    587                      &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   & 
    588                      &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
    589  
    590                   za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + & 
    591                      (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj)) 
    592                   zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + & 
    593                      za*(u_oce1(ji,jj)) 
    594                   zcca         = z0+za 
     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 
    595531                  zccb         = zcorl1(ji,jj) 
    596                   u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+epsd)*zmask  
    597                END DO ! ji 
    598             END DO ! jj 
     532                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
     533               END DO 
     534            END DO 
    599535 
    600536            CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    601537#if defined key_agrif && defined key_lim2 
    602             CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
     538            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
    603539#endif 
    604540#if defined key_bdy 
     
    611547            !---  Convergence test. 
    612548            DO jj = k_j1+1 , k_jpj-1 
    613                zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           & 
    614                   ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    615             END DO 
    616             zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) ) 
     549               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
     550            END DO 
     551            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) ) 
    617552            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
    618553         ENDIF 
     
    625560      ! 4) Prevent ice velocities when the ice is thin 
    626561      !------------------------------------------------------------------------------! 
    627       ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 
    628       ! ocean velocity,  
    629       ! This prevents high velocity when ice is thin 
    630 !CDIR NOVERRCHK 
     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 
    631564      DO jj = k_j1+1, k_jpj-1 
    632 !CDIR NOVERRCHK 
    633565         DO ji = fs_2, fs_jpim1 
    634             zdummy = vt_i(ji,jj) 
    635             IF ( zdummy .LE. hminrhg ) THEN 
     566            IF ( vt_i(ji,jj) <= zvmin ) THEN 
    636567               u_ice(ji,jj) = u_oce(ji,jj) 
    637568               v_ice(ji,jj) = v_oce(ji,jj) 
    638             ENDIF ! zdummy 
     569            ENDIF 
    639570         END DO 
    640571      END DO 
    641572 
    642       CALL lbc_lnk( u_ice(:,:), 'U', -1. )  
    643       CALL lbc_lnk( v_ice(:,:), 'V', -1. )  
     573      CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 
     574 
    644575#if defined key_agrif && defined key_lim2 
    645       CALL agrif_rhg_lim2( nevp , nevp, 'U' ) 
    646       CALL agrif_rhg_lim2( nevp , nevp, 'V' ) 
     576      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 
     577      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 
    647578#endif 
    648579#if defined key_bdy 
     
    653584      DO jj = k_j1+1, k_jpj-1  
    654585         DO ji = fs_2, fs_jpim1 
    655             zdummy = vt_i(ji,jj) 
    656             IF ( zdummy .LE. hminrhg ) THEN 
    657                v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   & 
    658                   &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 
    659                   &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
    660  
    661                u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   & 
    662                   &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 
    663                   &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
    664             ENDIF ! zdummy 
     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  
    665595         END DO 
    666596      END DO 
    667597 
    668       CALL lbc_lnk( u_ice2(:,:), 'V', -1. )  
    669       CALL lbc_lnk( v_ice1(:,:), 'U', -1. ) 
     598      CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 
    670599 
    671600      ! Recompute delta, shear and div, inputs for mechanical redistribution  
    672 !CDIR NOVERRCHK 
    673601      DO jj = k_j1+1, k_jpj-1 
    674 !CDIR NOVERRCHK 
    675          DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi 
     602         DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask 
    676603            !- divu_i(:,:), zdt(:,:): divergence and tension at centre  
    677604            !- zds(:,:): shear on northeast corner of grid cells 
    678             zdummy = vt_i(ji,jj) 
    679             IF ( zdummy .LE. hminrhg ) THEN 
    680  
    681                divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
    682                   &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  & 
    683                   &             +e1v(ji,jj)*v_ice(ji,jj)                      & 
    684                   &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  & 
    685                   &            )                                              & 
    686                   &            / area(ji,jj) 
    687  
    688                zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    & 
    689                   &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                & 
    690                   &           )*e2t(ji,jj)*e2t(ji,jj)                      & 
    691                   &          -( v_ice(ji,jj)/e1v(ji,jj)                    & 
    692                   &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                & 
    693                   &           )*e1t(ji,jj)*e1t(ji,jj)                      & 
    694                   &         )                                              & 
    695                   &        / area(ji,jj) 
     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_e12t(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_e12t(ji,jj) 
    696614               ! 
    697615               ! SB modif because ocean has no slip boundary condition  
    698                zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1)              & 
    699                   &           - u_ice(ji,jj)   / e1u(ji,jj) )              & 
    700                   &           * e1f(ji,jj) * e1f(ji,jj)                    & 
    701                   &          + ( v_ice(ji+1,jj) / e2v(ji+1,jj)             & 
    702                   &            - v_ice(ji,jj)  / e2v(ji,jj) )              & 
    703                   &           * e2f(ji,jj) * e2f(ji,jj) )                  & 
    704                   &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 
    705                   &        * tmi(ji,jj) * tmi(ji,jj+1)                     & 
    706                   &        * tmi(ji+1,jj) * tmi(ji+1,jj+1) 
    707  
    708                zdst = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    & 
    709                   &           - e2u( ji-1, jj   ) * v_ice1(ji-1,jj  )    & 
    710                   &           + e1v( ji  , jj   ) * u_ice2(ji  ,jj  )    & 
    711                   &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj) 
    712  
    713                delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )   
    714                delta_i(ji,jj) = delta + creepl 
     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_e12f(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_e12t(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 
    715626             
    716             ENDIF ! zdummy 
    717  
    718          END DO !jj 
    719       END DO !ji 
     627            ENDIF 
     628         END DO 
     629      END DO 
    720630      ! 
    721631      !------------------------------------------------------------------------------! 
    722632      ! 5) Store stress tensor and its invariants 
    723633      !------------------------------------------------------------------------------! 
    724       ! 
    725634      ! * Invariants of the stress tensor are required for limitd_me 
    726635      !   (accelerates convergence and improves stability) 
    727636      DO jj = k_j1+1, k_jpj-1 
    728637         DO ji = fs_2, fs_jpim1 
    729             ! begin TECLIM change  
    730             zdst= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &    
    731                &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &    
    732                &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &    
    733                &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj)  
     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_e12t(ji,jj)  
    734640            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 
    735             ! end TECLIM change 
    736641         END DO 
    737642      END DO 
    738643 
    739644      ! Lateral boundary condition 
    740       CALL lbc_lnk( divu_i (:,:), 'T', 1. ) 
    741       CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 
    742       ! CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 
    743       CALL lbc_lnk( shear_i(:,:), 'T', 1. ) 
     645      CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. ) 
    744646 
    745647      ! * Store the stress tensor for the next time step 
     
    772674            DO jj = k_j1+1, k_jpj-1 
    773675               DO ji = 2, jpim1 
    774                   IF (zpresh(ji,jj) .GT. 1.0) THEN 
     676                  IF (zpresh(ji,jj) > 1.0) THEN 
    775677                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )  
    776678                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
     
    786688      ! 
    787689      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    788       CALL wrk_dealloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                ) 
     690      CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    789691      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    790692      CALL wrk_dealloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
Note: See TracChangeset for help on using the changeset viewer.