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 11692 for NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/ICE/icedyn_rhg_evp.F90 – NEMO

Ignore:
Timestamp:
2019-10-12T16:08:18+02:00 (4 years ago)
Author:
francesca
Message:

Update branch to integrate the development starting from the current v4.01 ready trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/src/ICE/icedyn_rhg_evp.F90

    r10891 r11692  
    112112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    113113      !! 
    114       LOGICAL, PARAMETER ::   ll_bdy_substep = .TRUE. ! temporary option to call bdy at each sub-time step (T) 
    115       !                                                                              or only at the main time step (F) 
    116114      INTEGER ::   ji, jj       ! dummy loop indices 
    117115      INTEGER ::   jter         ! local integers 
     
    123121      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    124122      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars 
    125       REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars 
     123      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars 
    126124      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    127125      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast 
     
    132130      REAL(wp) ::   zshear, zdum1, zdum2 
    133131      ! 
    134       REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors 
    135132      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    136133      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    137134      ! 
    138135      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points 
    139       REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points 
     136      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points 
    140137      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    141138      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    142       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
    143       REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ib , ztauV_ib             ! ice-bottom stress at U-V points (landfast param) 
    144       REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points 
    145139      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
    146       REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses 
    147140      ! 
    148141      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
     
    152145      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
    153146      !                                                                 !    ice bottom surface if ice is embedded    
    154       REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array 
    155       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array 
    156       ! 
    157       REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays 
    158       REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence 
     147      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses 
     148      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points 
     149      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array 
     150      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points 
     151      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points 
     152      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast) 
     153      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast) 
     154      ! 
     155      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays 
     156      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence 
    159157      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice 
    160158 
     
    163161      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    164162      !! --- diags 
    165       REAL(wp), DIMENSION(jpi,jpj) ::   zswi 
     163      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00 
    166164      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3 
    167165      !! --- SIMIP diags 
    168       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice    
    169       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice 
    170       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2) 
    171       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2) 
    172       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2) 
    173       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2) 
    174       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2) 
    175       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2) 
    176       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress 
    177       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress   
    178166      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
    179167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) 
     
    255243      CALL ice_strength 
    256244 
    257       ! scale factors 
    258       DO jj = 2, jpjm1 
    259          DO ji = fs_2, fs_jpim1 
    260             z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) ) 
    261             z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) ) 
    262          END DO 
    263       END DO 
    264  
    265245      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
    266       IF( ln_landfast_L16 .OR. ln_landfast_home ) THEN   ;   zkt = rn_tensile 
    267       ELSE                                               ;   zkt = 0._wp 
     246      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile 
     247      ELSE                         ;   zkt = 0._wp 
    268248      ENDIF 
    269249      ! 
     
    291271 
    292272            ! Ocean currents at U-V points 
    293             v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    & 
    294                &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    295              
    296             u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    & 
    297                &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     273            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 
     274            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 
    298275 
    299276            ! Coriolis at T points (m*f) 
     
    308285             
    309286            ! Drag ice-atm. 
    310             zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
    311             zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
     287            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     288            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
    312289 
    313290            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
     
    316293 
    317294            ! masks 
    318             zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
    319             zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
     295            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
     296            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
    320297 
    321298            ! switches 
    322             IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zswitchU(ji,jj) = 0._wp 
    323             ELSE                                                   ;   zswitchU(ji,jj) = 1._wp   ;   ENDIF 
    324             IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zswitchV(ji,jj) = 0._wp 
    325             ELSE                                                   ;   zswitchV(ji,jj) = 1._wp   ;   ENDIF 
     299            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
     300            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
     301            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
     302            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF 
    326303 
    327304         END DO 
     
    339316               ! ice-bottom stress at U points 
    340317               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 
    341                zTauU_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
     318               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    342319               ! ice-bottom stress at V points 
    343320               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 
    344                zTauV_ib(ji,jj)   = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
     321               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    345322               ! ice_bottom stress at T points 
    346323               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 
    347                tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
     324               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    348325            END DO 
    349326         END DO 
    350327         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1. ) 
    351328         ! 
    352       ELSEIF( ln_landfast_home ) THEN          !-- Home made 
     329      ELSE                               !-- no landfast 
    353330         DO jj = 2, jpjm1 
    354331            DO ji = fs_2, fs_jpim1 
    355                zTauU_ib(ji,jj) = tau_icebfr(ji,jj) 
    356                zTauV_ib(ji,jj) = tau_icebfr(ji,jj) 
    357             END DO 
    358          END DO 
    359          ! 
    360       ELSE                                     !-- no landfast 
    361          DO jj = 2, jpjm1 
    362             DO ji = fs_2, fs_jpim1 
    363                zTauU_ib(ji,jj) = 0._wp 
    364                zTauV_ib(ji,jj) = 0._wp 
     332               ztaux_base(ji,jj) = 0._wp 
     333               ztauy_base(ji,jj) = 0._wp 
    365334            END DO 
    366335         END DO 
    367336      ENDIF 
    368       IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr(:,:) ) 
    369337 
    370338      !------------------------------------------------------------------------------! 
     
    372340      !------------------------------------------------------------------------------! 
    373341      ! 
    374       !                                               !----------------------! 
     342      !                                               ! ==================== ! 
    375343      DO jter = 1 , nn_nevp                           !    loop over jter    ! 
    376          !                                            !----------------------!         
     344         !                                            ! ==================== !         
    377345         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
    378346         ! 
     
    479447                  &                  ) * r1_e1e2v(ji,jj) 
    480448               ! 
    481                !                !--- u_ice at V point 
    482                u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    483                   &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     449               !                !--- ice currents at U-V point 
     450               v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1) 
     451               u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1) 
    484452               ! 
    485                !                !--- v_ice at U point 
    486                v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    487                   &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    488453            END DO 
    489454         END DO 
     
    504469                  !                 !--- tau_bottom/v_ice 
    505470                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    506                   zTauB = - zTauV_ib(ji,jj) / zvel 
     471                  zTauB = ztauy_base(ji,jj) / zvel 
     472                  !                 !--- OceanBottom-to-Ice stress 
     473                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    507474                  ! 
    508475                  !                 !--- Coriolis at V-points (energy conserving formulation) 
    509                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     476                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    510477                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    511478                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    512479                  ! 
    513480                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    514                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    515                   ! 
    516                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    517                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     481                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     482                  ! 
     483                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     484                  !                                         1 = sliding friction : TauB < RHS 
     485                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    518486                  ! 
    519487                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    520                   v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
    521                      &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    522                      &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    523                                     + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    524                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    525                      &           ) * zmaskV(ji,jj) 
     488                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
     489                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     490                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     491                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     492                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     493                        &           )   * zmsk00y(ji,jj) 
    526494                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    527                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             & ! previous velocity 
    528                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    529                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    530                      &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    531                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    532                      &            ) * zmaskV(ji,jj) 
     495                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
     496                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     497                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     498                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     499                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     500                        &            )   * zmsk00y(ji,jj) 
    533501                  ENDIF 
    534502               END DO 
     
    540508            CALL agrif_interp_ice( 'V' ) 
    541509#endif 
    542             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' ) 
     510            IF( ln_bdy CALL bdy_ice_dyn( 'V' ) 
    543511            ! 
    544512            DO jj = 2, jpjm1 
     
    552520                  !                 !--- tau_bottom/u_ice 
    553521                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    554                   zTauB = - zTauU_ib(ji,jj) / zvel 
     522                  zTauB = ztaux_base(ji,jj) / zvel 
     523                  !                 !--- OceanBottom-to-Ice stress 
     524                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
    555525                  ! 
    556526                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    557                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
     527                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    558528                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    559529                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    560530                  ! 
    561531                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    562                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    563                   ! 
    564                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    565                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     532                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     533                  ! 
     534                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     535                  !                                         1 = sliding friction : TauB < RHS 
     536                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    566537                  ! 
    567538                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    568                   u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
    569                      &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    570                      &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    571                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    572                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    573                      &            ) * zmaskU(ji,jj) 
     539                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
     540                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     541                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     542                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     543                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     544                        &           )   * zmsk00x(ji,jj) 
    574545                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    575                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             & ! previous velocity 
    576                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    577                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    578                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    579                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    580                      &            ) * zmaskU(ji,jj) 
     546                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
     547                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     548                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     549                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     550                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     551                        &            )   * zmsk00x(ji,jj) 
    581552                  ENDIF 
    582553               END DO 
     
    588559            CALL agrif_interp_ice( 'U' ) 
    589560#endif 
    590             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' ) 
     561            IF( ln_bdy CALL bdy_ice_dyn( 'U' ) 
    591562            ! 
    592563         ELSE ! odd iterations 
     
    602573                  !                 !--- tau_bottom/u_ice 
    603574                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) 
    604                   zTauB = - zTauU_ib(ji,jj) / zvel 
     575                  zTauB = ztaux_base(ji,jj) / zvel 
     576                  !                 !--- OceanBottom-to-Ice stress 
     577                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj) 
    605578                  ! 
    606579                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    607                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
     580                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    608581                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    609582                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    610583                  ! 
    611584                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    612                   zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    613                   ! 
    614                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    615                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - zTauU_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     585                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
     586                  ! 
     587                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     588                  !                                         1 = sliding friction : TauB < RHS 
     589                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    616590                  ! 
    617591                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    618                   u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
    619                      &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    620                      &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    621                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
    622                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    623                      &            ) * zmaskU(ji,jj) 
     592                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
     593                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     594                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     595                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     596                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     597                        &           )   * zmsk00x(ji,jj) 
    624598                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    625                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             & ! previous velocity 
    626                      &                                     + zTauE + zTauO * u_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    627                      &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    628                      &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    629                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchU(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    630                      &            ) * zmaskU(ji,jj) 
     599                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity 
     600                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     601                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     602                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     603                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     604                        &            )   * zmsk00x(ji,jj) 
    631605                  ENDIF 
    632606               END DO 
     
    638612            CALL agrif_interp_ice( 'U' ) 
    639613#endif 
    640             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'U' ) 
     614            IF( ln_bdy CALL bdy_ice_dyn( 'U' ) 
    641615            ! 
    642616            DO jj = 2, jpjm1 
     
    650624                  !                 !--- tau_bottom/v_ice 
    651625                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 
    652                   zTauB = - zTauV_ib(ji,jj) / zvel 
     626                  zTauB = ztauy_base(ji,jj) / zvel 
     627                  !                 !--- OceanBottom-to-Ice stress 
     628                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 
    653629                  ! 
    654630                  !                 !--- Coriolis at v-points (energy conserving formulation) 
    655                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     631                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    656632                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    657633                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    658634                  ! 
    659635                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    660                   zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    661                   ! 
    662                   !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    663                   rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - zTauV_ib(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
     636                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
     637                  ! 
     638                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS) 
     639                  !                                         1 = sliding friction : TauB < RHS 
     640                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) ) 
    664641                  ! 
    665642                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    666                   v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
    667                      &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    668                      &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    669                                     + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    670                      &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    671                      &           ) * zmaskV(ji,jj) 
     643                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
     644                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     645                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     646                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     647                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     648                        &           )   * zmsk00y(ji,jj) 
    672649                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    673                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             & ! previous velocity 
    674                      &                                     + zTauE + zTauO * v_ice(ji,jj)                             & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    675                      &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             & ! m/dt + tau_io(only ice part) + landfast 
    676                      &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  & ! static friction => slow decrease to v=0 
    677                      &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zswitchV(ji,jj) )        & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    678                      &            ) * zmaskV(ji,jj) 
     650                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity 
     651                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     652                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
     653                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
     654                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     655                        &            )   * zmsk00y(ji,jj) 
    679656                  ENDIF 
    680657               END DO 
     
    686663            CALL agrif_interp_ice( 'V' ) 
    687664#endif 
    688             IF( ln_bdy .AND. ll_bdy_substep ) CALL bdy_ice_dyn( 'V' ) 
     665            IF( ln_bdy CALL bdy_ice_dyn( 'V' ) 
    689666            ! 
    690667         ENDIF 
     
    701678      END DO                                              !  end loop over jter  ! 
    702679      !                                                   ! ==================== ! 
    703       ! 
    704       IF( ln_bdy .AND. .NOT.ll_bdy_substep ) THEN 
    705          CALL bdy_ice_dyn( 'U' ) 
    706          CALL bdy_ice_dyn( 'V' ) 
    707       ENDIF 
    708680      ! 
    709681      !------------------------------------------------------------------------------! 
     
    764736      DO jj = 1, jpj 
    765737         DO ji = 1, jpi 
    766             zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
     738            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    767739         END DO 
    768740      END DO 
    769741 
     742      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
     743      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 
     744         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 
     745         ! 
     746         CALL lbc_lnk_multi( 'icedyn_rhg_evp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., & 
     747            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 
     748         ! 
     749         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 
     750         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 ) 
     751         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 ) 
     752         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 ) 
     753         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 ) 
     754         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 ) 
     755      ENDIF 
     756        
    770757      ! --- divergence, shear and strength --- ! 
    771       IF( iom_use('icediv') )   CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) )   ! divergence 
    772       IF( iom_use('iceshe') )   CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:) )   ! shear 
    773       IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength 
    774  
    775       ! --- charge ellipse --- ! 
    776       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN 
     758      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
     759      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
     760      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
     761 
     762      ! --- stress tensor --- ! 
     763      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    777764         ! 
    778765         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     
    780767         DO jj = 2, jpjm1 
    781768            DO ji = 2, jpim1 
    782                zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    783                   &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
    784                   &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 
     769               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
     770                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
     771                  &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
    785772 
    786773               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    787774 
    788                zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
     775               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    789776 
    790777!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
     
    799786         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
    800787         ! 
    801          IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 ) 
    802          IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 ) 
    803          IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 ) 
    804          ! 
     788         CALL iom_put( 'isig1' , zsig1 ) 
     789         CALL iom_put( 'isig2' , zsig2 ) 
     790         CALL iom_put( 'isig3' , zsig3 ) 
     791         ! 
     792         ! Stress tensor invariants (normal and shear stress N/m) 
     793         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
     794         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
     795 
    805796         DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
    806797      ENDIF 
    807798       
    808799      ! --- SIMIP --- ! 
    809       IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. & 
    810          & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. & 
    811          & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & 
    812          & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN 
    813  
    814          ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  & 
    815             &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  & 
    816             &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  & 
    817             &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) ) 
    818           
     800      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
     801         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
     802         ! 
     803         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zspgU, 'U', -1., zspgV, 'V', -1., & 
     804            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 
     805 
     806         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x) 
     807         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y) 
     808         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x) 
     809         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y) 
     810         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x) 
     811         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y) 
     812      ENDIF 
     813 
     814      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. & 
     815         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN 
     816         ! 
     817         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , & 
     818            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 
     819         ! 
    819820         DO jj = 2, jpjm1 
    820821            DO ji = 2, jpim1 
    821                rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    822                 
    823                ! Stress tensor invariants (normal and shear stress N/m) 
    824                zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress 
    825                zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress 
    826                 
    827                ! Stress terms of the momentum equation (N/m2) 
    828                zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch     ! sea surface slope stress term 
    829                zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 
    830                 
    831                zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch     ! Coriolis stress term 
    832                zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch 
    833                 
    834                zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch     ! internal stress term 
    835                zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch 
    836                 
    837                zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress 
    838                zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 
    839                 
    840822               ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    841                zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 
    842                zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
    843                 
     823               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 
     824               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 
     825 
    844826               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
    845827               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
    846                 
     828 
    847829               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
    848830               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
    849                 
     831 
    850832               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component 
    851833               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   '' 
    852                 
    853             END DO 
    854          END DO 
    855           
    856          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   & 
    857             &                zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   & 
    858             &                zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   &  
    859             &                zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    ) 
    860                    
    861          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi  , 'U', -1., zdiag_vtau_oi  , 'V', -1.,   & 
    862             &                zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   & 
    863             &                zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   & 
    864             &                zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    ) 
    865           
    866          IF( iom_use('normstr' ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress 
    867          IF( iom_use('sheastr' ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress 
    868          IF( iom_use('dssh_dx' ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x) 
    869          IF( iom_use('dssh_dy' ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y) 
    870          IF( iom_use('corstrx' ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x) 
    871          IF( iom_use('corstry' ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y) 
    872          IF( iom_use('intstrx' ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x) 
    873          IF( iom_use('intstry' ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y) 
    874          IF( iom_use('utau_oi' ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x) 
    875          IF( iom_use('vtau_oi' ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y) 
    876          IF( iom_use('xmtrpice') )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s) 
    877          IF( iom_use('ymtrpice') )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport  
    878          IF( iom_use('xmtrpsnw') )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s) 
    879          IF( iom_use('ymtrpsnw') )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport 
    880          IF( iom_use('xatrp'   ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport 
    881          IF( iom_use('yatrp'   ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport 
    882  
    883          DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  & 
    884             &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  & 
    885             &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  & 
    886             &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     ) 
     834 
     835            END DO 
     836         END DO 
     837 
     838         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 
     839            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & 
     840            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. ) 
     841 
     842         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s) 
     843         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport  
     844         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s) 
     845         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport 
     846         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport 
     847         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport 
     848 
     849         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , & 
     850            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp ) 
    887851 
    888852      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.