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 13746 for NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90 – NEMO

Ignore:
Timestamp:
2020-11-08T22:25:28+01:00 (3 years ago)
Author:
stefryn
Message:

update test case

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_evp.F90

    r13574 r13746  
    4141   USE prtctl         ! Print control 
    4242 
     43   USE netcdf         ! NetCDF library for convergence test 
    4344   IMPLICIT NONE 
    4445   PRIVATE 
     
    4950   !! * Substitutions 
    5051#  include "vectopt_loop_substitute.h90" 
     52 
     53   !! for convergence tests 
     54   INTEGER ::   ncvgid   ! netcdf file id 
     55   INTEGER ::   nvarid   ! netcdf variable id 
     56   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
    5157   !!---------------------------------------------------------------------- 
    5258   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    53    !! $Id: icedyn_rhg_evp.F90 11536 2019-09-11 13:54:18Z smasson $ 
     59   !! $Id: icedyn_rhg_evp.F90 13662 2020-10-22 18:49:56Z clem $ 
    5460   !! Software governed by the CeCILL license (see ./LICENSE) 
    5561   !!---------------------------------------------------------------------- 
     
    119125      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity 
    120126      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017 
     127      REAl(wp) ::   zbetau, zbetav 
    121128      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    122       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars 
     129      REAL(wp) ::   zp_delf, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars 
    123130      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars 
    124131      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    125132      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast 
    126133      ! 
    127       REAL(wp) ::   zresm                                               ! Maximal error on ice velocity 
    128134      REAL(wp) ::   zintb, zintn                                        ! dummy argument 
    129135      REAL(wp) ::   zfac_x, zfac_y 
    130136      REAL(wp) ::   zshear, zdum1, zdum2 
    131       REAL(wp) ::   invw                                                ! for test case 
    132       ! 
    133       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
     137      REAL(wp) ::   zinvw                                                ! for test case 
     138 
     139      ! 
     140      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points 
     141      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension 
    134142      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    135143      ! 
     
    138146      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    139147      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    140       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
     148      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points 
    141149      ! 
    142150      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
    143151      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
    144 !!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence 
    145152      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
    146153      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
     
    156163      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays 
    157164      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence 
    158       REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice 
     165      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice 
    159166 
    160167      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
    161168      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small 
    162169      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
     170      !! --- check convergence 
     171      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice 
    163172      !! --- diags 
    164       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00 
    165       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3 
     173      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength 
     174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p          
    166175      !! --- SIMIP diags 
    167176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
     
    175184      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' 
    176185      ! 
    177 !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization.... 
     186      ! for diagnostics and convergence tests 
     187      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 
     188      DO jj = 1, jpj 
     189         DO ji = 1, jpi 
     190            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice 
     191            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 
     192         END DO 
     193      END DO 
     194      ! 
     195      !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization.... 
    178196      !------------------------------------------------------------------------------! 
    179197      ! 0) mask at F points for the ice 
     
    188206 
    189207      ! Lateral boundary conditions on velocity (modify zfmask) 
    190       zwf(:,:) = zfmask(:,:) 
    191208      DO jj = 2, jpjm1 
    192209         DO ji = fs_2, fs_jpim1   ! vector opt. 
    193210            IF( zfmask(ji,jj) == 0._wp ) THEN 
    194                zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 
     211               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 
     212                  &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 
    195213            ENDIF 
    196214         END DO 
     
    198216      DO jj = 2, jpjm1 
    199217         IF( zfmask(1,jj) == 0._wp ) THEN 
    200             zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 
     218            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 
    201219         ENDIF 
    202220         IF( zfmask(jpi,jj) == 0._wp ) THEN 
    203             zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
    204          ENDIF 
     221            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 
     222        ENDIF 
    205223      END DO 
    206224      DO ji = 2, jpim1 
    207225         IF( zfmask(ji,1) == 0._wp ) THEN 
    208             zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 
     226            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 
    209227         ENDIF 
    210228         IF( zfmask(ji,jpj) == 0._wp ) THEN 
    211             zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 
     229            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 
    212230         ENDIF 
    213231      END DO 
     
    219237      zrhoco = rau0 * rn_cio  
    220238!extra code for test case boundary conditions 
    221       invw=1._wp/(zrhoco*0.5_wp) 
     239      zinvw=1._wp/(zrhoco*0.5_wp) 
    222240 
    223241      ! ecc2: square of yield ellipse eccenticrity 
     
    225243      z1_ecc2 = 1._wp / ecc2 
    226244 
    227       ! Time step for subcycling 
    228       zdtevp   = rdt_ice / REAL( nn_nevp ) 
    229       z1_dtevp = 1._wp / zdtevp 
    230  
    231245      ! alpha parameters (Bouillon 2009) 
    232246      IF( .NOT. ln_aEVP ) THEN 
    233          zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
     247         zdtevp   = rdt_ice / REAL( nn_nevp ) 
     248         zalph1 =   2._wp * rn_relast * REAL( nn_nevp ) 
    234249         zalph2 = zalph1 * z1_ecc2 
    235250 
    236251         z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    237252         z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     253      ELSE 
     254         zdtevp   = rdt_ice 
     255         ! zalpha parameters set later on adaptatively 
    238256      ENDIF 
     257      z1_dtevp = 1._wp / zdtevp 
    239258          
    240259      ! Initialise stress tensor  
     
    247266 
    248267      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
    249       IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile 
     268      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile 
    250269      ELSE                         ;   zkt = 0._wp 
    251270      ENDIF 
     
    318337               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    319338               ! ice-bottom stress at U points 
    320                zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj) 
    321                ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
     339               zvCr = zaU(ji,jj) * rn_lf_depfra * hu_n(ji,jj) 
     340               ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    322341               ! ice-bottom stress at V points 
    323                zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj) 
    324                ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
     342               zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj) 
     343               ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    325344               ! ice_bottom stress at T points 
    326                zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj) 
    327                tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
     345               zvCr = at_i(ji,jj) * rn_lf_depfra * ht_n(ji,jj) 
     346               tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    328347            END DO 
    329348         END DO 
     
    348367         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
    349368         ! 
    350 !!$         IF(ln_ctl) THEN   ! Convergence test 
    351 !!$            DO jj = 1, jpjm1 
    352 !!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
    353 !!$               zv_ice(:,jj) = v_ice(:,jj) 
    354 !!$            END DO 
    355 !!$         ENDIF 
     369         ! convergence test 
     370         IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2  ) THEN 
     371            DO jj = 1, jpj 
     372               DO ji = 1, jpi 
     373                  zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 
     374                  zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 
     375               END DO 
     376            END DO 
     377         ENDIF 
    356378 
    357379         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
    358          DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
     380         DO jj = 1, jpjm1 
    359381            DO ji = 1, jpim1 
    360382 
     
    366388            END DO 
    367389         END DO 
    368          CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. ) 
    369  
    370          DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
    371             DO ji = 2, jpi ! no vector loop 
     390 
     391         DO jj = 2, jpjm1 
     392            DO ji = 2, jpim1 ! no vector loop 
    372393 
    373394               ! shear**2 at T points (doc eq. A16) 
     
    389410                
    390411               ! delta at T points 
    391                zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
    392  
    393                ! P/delta at T points 
    394                zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
    395  
    396                ! alpha & beta for aEVP 
     412               zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     413 
     414            END DO 
     415         END DO 
     416         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1._wp ) 
     417          
     418         ! P/delta at T points 
     419         DO jj = 1, jpj 
     420            DO ji = 1, jpi 
     421               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 
     422            END DO 
     423         END DO 
     424 
     425         DO jj = 2, jpj    ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 
     426            DO ji = 2, jpi ! no vector loop 
     427 
     428               ! divergence at T points (duplication to avoid communications) 
     429               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     430                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     431                  &    ) * r1_e1e2t(ji,jj) 
     432                
     433               ! tension at T points (duplication to avoid communications) 
     434               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     435                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     436                  &   ) * r1_e1e2t(ji,jj) 
     437 
     438               ! alpha for aEVP 
    397439               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
    398440               !   alpha = beta = sqrt(4*gamma) 
     
    402444                  zalph2   = zalph1 
    403445                  z1_alph2 = z1_alph1 
     446                  ! explicit: 
     447                  ! z1_alph1 = 1._wp / zalph1 
     448                  ! z1_alph2 = 1._wp / zalph1 
     449                  ! zalph1 = zalph1 - 1._wp 
     450                  ! zalph2 = zalph1 
    404451               ENDIF 
    405452                
    406453               ! stress at T points (zkt/=0 if landfast) 
    407                zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
    408                zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
     454               zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 
     455               zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
    409456              
    410457            END DO 
    411458         END DO 
    412          CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. ) 
    413  
     459 
     460         ! Save beta at T-points for further computations 
     461         IF( ln_aEVP ) THEN 
     462            DO jj = 1, jpj 
     463               DO ji = 1, jpi 
     464                  zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     465               END DO 
     466            END DO 
     467         ENDIF 
     468          
    414469         DO jj = 1, jpjm1 
    415470            DO ji = 1, jpim1 
    416471 
    417                ! alpha & beta for aEVP 
     472               ! alpha for aEVP 
    418473               IF( ln_aEVP ) THEN 
    419                   zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     474                  zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 
    420475                  z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
    421                   zbeta(ji,jj) = zalph2 
     476                  ! explicit: 
     477                  ! z1_alph2 = 1._wp / zalph2 
     478                  ! zalph2 = zalph2 - 1._wp 
    422479               ENDIF 
    423480                
     
    489546                  ! 
    490547                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    491                      v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    492                         &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    493                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    494                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    495                         &             ) * 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 
     548                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 
     549                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     550                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     551                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     552                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &  
     553                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     554                        &                                    ) / ( zbetav + 1._wp )                                              & 
     555                        &             ) * 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 
    496556                        &           )   * zmsk00y(ji,jj) 
    497557                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    498                      v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
    499                         &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    500                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    501                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    502                         &              ) * 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 
    503                         &            )   * zmsk00y(ji,jj) 
     558                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
     559                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     560                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     561                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     562                        &             ) * 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 
     563                        &            )  * zmsk00y(ji,jj) 
    504564                  ENDIF 
    505565!extra code for test case boundary conditions 
    506566                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    507                      v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 
     567                     v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 
    508568                  END IF 
    509569               END DO 
     
    544604                  ! 
    545605                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    546                      u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(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) * ( zbeta(ji,jj) + 1._wp ) + 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  
     606                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 
     607                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     608                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     609                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     610                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   & 
     611                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     612                        &                                    ) / ( zbetau + 1._wp )                                              & 
     613                        &             ) * 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  
    551614                        &           )   * zmsk00x(ji,jj) 
    552615                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    553                      u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
    554                         &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    555                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    556                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    557                         &              ) * 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  
    558                         &            )   * zmsk00x(ji,jj) 
     616                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
     617                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     618                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     619                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     620                        &             ) * 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  
     621                        &           )   * zmsk00x(ji,jj) 
    559622                  ENDIF 
    560623!extra code for test case boundary conditions 
    561624                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    562                      u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 
     625                     u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 
    563626                  END IF 
    564627               END DO 
     
    601664                  ! 
    602665                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    603                      u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    604                         &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    605                         &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    606                         &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    607                         &             ) * 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  
     666                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 
     667                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     668                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     669                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     670                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   & 
     671                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     672                        &                                    ) / ( zbetau + 1._wp )                                              & 
     673                        &             ) * 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  
    608674                        &           )   * zmsk00x(ji,jj) 
    609675                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    610                      u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
    611                         &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    612                         &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    613                         &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    614                         &              ) * 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 
    615                         &            )   * zmsk00x(ji,jj) 
     676                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
     677                        &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     678                        &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     679                        &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     680                        &             ) * 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 
     681                        &           )   * zmsk00x(ji,jj) 
    616682                  ENDIF 
    617683!extra code for test case boundary conditions 
    618684                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    619                      u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 
     685                     u_ice(ji,jj) = zinvw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 
    620686                  END IF 
    621687               END DO 
     
    656722                  ! 
    657723                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    658                      v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    659                         &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    660                         &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    661                         &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    662                         &             ) * 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 
     724                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 
     725                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     726                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     727                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     728                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
     729                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     730                        &                                    ) / ( zbetav + 1._wp )                                              &  
     731                        &             ) * 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 
    663732                        &           )   * zmsk00y(ji,jj) 
    664733                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    665                      v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
    666                         &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    667                         &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    668                         &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    669                         &              ) * 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 
    670                         &            )   * zmsk00y(ji,jj) 
     734                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
     735                        &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     736                        &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     737                        &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     738                        &             ) * 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 
     739                        &           )   * zmsk00y(ji,jj) 
    671740                  ENDIF 
    672741!extra code for test case boundary conditions 
    673742                  IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 
    674                      v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 
     743                     v_ice(ji,jj) = zinvw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 
    675744                  END IF 
    676745               END DO 
     
    686755         ENDIF 
    687756 
    688 !!$         IF(ln_ctl) THEN   ! Convergence test 
    689 !!$            DO jj = 2 , jpjm1 
    690 !!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    691 !!$            END DO 
    692 !!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 
    693 !!$            CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
    694 !!$         ENDIF 
     757         ! convergence test 
     758         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 
    695759         ! 
    696760         !                                                ! ==================== ! 
    697761      END DO                                              !  end loop over jter  ! 
    698762      !                                                   ! ==================== ! 
     763      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta ) 
    699764      ! 
    700765      !------------------------------------------------------------------------------! 
     
    721786            zdt2 = zdt * zdt 
    722787             
     788            zten_i(ji,jj) = zdt 
     789 
    723790            ! shear**2 at T points (doc eq. A16) 
    724791            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     
    735802             
    736803            ! delta at T points 
    737             zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    738             rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    739             pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     804            zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta   
     805            rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0 
     806            pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl 
    740807 
    741808         END DO 
    742809      END DO 
    743       CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
     810      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., & 
     811         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1. ) 
    744812       
    745813      ! --- Store the stress tensor for the next time step --- ! 
    746       CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    747814      pstress1_i (:,:) = zs1 (:,:) 
    748815      pstress2_i (:,:) = zs2 (:,:) 
     
    753820      ! 5) diagnostics 
    754821      !------------------------------------------------------------------------------! 
    755       DO jj = 1, jpj 
    756          DO ji = 1, jpi 
    757             zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    758          END DO 
    759       END DO 
    760  
    761822      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
    762823      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 
     
    777838      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
    778839      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
     840      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta 
    779841      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    780842 
    781       ! --- stress tensor --- ! 
    782       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    783          ! 
    784          ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     843      ! --- Stress tensor invariants (SIMIP diags) --- ! 
     844      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
     845         ! 
     846         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    785847         !          
    786          DO jj = 2, jpjm1 
    787             DO ji = 2, jpim1 
    788                zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    789                   &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
    790                   &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) ) 
    791  
    792                zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
    793  
    794                zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    795  
    796 !!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    797 !!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
    798 !!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
    799 !!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    800                zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    801                zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
    802                zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    803             END DO 
    804          END DO 
    805          CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
    806          ! 
    807          CALL iom_put( 'isig1' , zsig1 ) 
    808          CALL iom_put( 'isig2' , zsig2 ) 
    809          CALL iom_put( 'isig3' , zsig3 ) 
    810          ! 
    811          ! Stress tensor invariants (normal and shear stress N/m) 
    812          IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress 
    813          IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
    814  
    815          DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
     848         DO jj = 1, jpj 
     849            DO ji = 1, jpi 
     850             
     851               ! Ice stresses 
     852               ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013) 
     853               ! These are NOT stress tensor components, neither stress invariants, neither stress principal components 
     854               ! I know, this can be confusing... 
     855               zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )  
     856               zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 
     857               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     858               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     859                
     860               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 
     861               zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
     862               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
     863                
     864            END DO 
     865         END DO          
     866         ! 
     867         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008) 
     868         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress 
     869         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress 
     870          
     871         DEALLOCATE ( zsig_I, zsig_II ) 
     872          
    816873      ENDIF 
     874 
     875      ! --- Normalized stress tensor principal components --- ! 
     876      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020 
     877      ! Recommendation 1 : we use ice strength, not replacement pressure 
     878      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities 
     879      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
     880         ! 
     881         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )          
     882         !          
     883         DO jj = 1, jpj 
     884            DO ji = 1, jpi 
     885             
     886               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
     887               !                        and **deformations** at current iterates 
     888               !                        following Lemieux & Dupont (2020) 
     889               zfac             =   zp_delt(ji,jj) 
     890               zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) ) 
     891               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
     892               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     893                
     894               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
     895               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure 
     896               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress 
    817897       
     898               ! Normalized  principal stresses (used to display the ellipse) 
     899               z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) ) 
     900               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
     901               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
     902            END DO 
     903         END DO                
     904         ! 
     905         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
     906         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
     907       
     908         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II ) 
     909          
     910      ENDIF 
     911 
    818912      ! --- SIMIP --- ! 
    819913      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
     
    871965      ENDIF 
    872966      ! 
     967      ! --- convergence tests --- ! 
     968      IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 
     969         IF( iom_use('uice_cvg') ) THEN 
     970            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
     971               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 
     972                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     973            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
     974               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 
     975                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     976            ENDIF 
     977         ENDIF 
     978      ENDIF       
     979      ! 
     980      DEALLOCATE( zmsk00, zmsk15 ) 
     981      ! 
    873982   END SUBROUTINE ice_dyn_rhg_evp 
     983 
     984 
     985   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 
     986      !!---------------------------------------------------------------------- 
     987      !!                    ***  ROUTINE rhg_cvg  *** 
     988      !!                      
     989      !! ** Purpose :   check convergence of oce rheology 
     990      !! 
     991      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity 
     992      !!                during the sub timestepping of rheology so as: 
     993      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 
     994      !!                This routine is called every sub-iteration, so it is cpu expensive 
     995      !! 
     996      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)    
     997      !!---------------------------------------------------------------------- 
     998      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index 
     999      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities 
     1000      !! 
     1001      INTEGER           ::   it, idtime, istatus 
     1002      INTEGER           ::   ji, jj          ! dummy loop indices 
     1003      REAL(wp)          ::   zresm           ! local real  
     1004      CHARACTER(len=20) ::   clname 
     1005      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence 
     1006      !!---------------------------------------------------------------------- 
     1007 
     1008      ! create file 
     1009      IF( kt == nit000 .AND. kiter == 1 ) THEN 
     1010         ! 
     1011         IF( lwp ) THEN 
     1012            WRITE(numout,*) 
     1013            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 
     1014            WRITE(numout,*) '~~~~~~~' 
     1015         ENDIF 
     1016         ! 
     1017         IF( lwm ) THEN 
     1018            clname = 'ice_cvg.nc' 
     1019            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 
     1020            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 
     1021            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime ) 
     1022            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid ) 
     1023            istatus = NF90_ENDDEF(ncvgid) 
     1024         ENDIF 
     1025         ! 
     1026      ENDIF 
     1027 
     1028      ! time 
     1029      it = ( kt - 1 ) * kitermax + kiter 
     1030       
     1031      ! convergence 
     1032      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 
     1033         zresm = 0._wp 
     1034      ELSE 
     1035         DO jj = 1, jpj 
     1036            DO ji = 1, jpi 
     1037               zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     1038                  &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 
     1039            END DO 
     1040         END DO 
     1041 
     1042         ! cut of the boundary of the box (forced velocities) 
     1043         IF (mjg(jj)<=30 .or. mjg(jj)>970 .or. mig(ji)<=30 .or. mig(ji)>970) THEN 
     1044            zres(ji,jj) = 0._wp 
     1045         END IF 
     1046 
     1047         zresm = MAXVAL( zres ) 
     1048         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
     1049      ENDIF 
     1050 
     1051      IF( lwm ) THEN 
     1052         ! write variables 
     1053         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 
     1054         ! close file 
     1055         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid) 
     1056      ENDIF 
     1057       
     1058   END SUBROUTINE rhg_cvg 
    8741059 
    8751060 
     
    9291114   END SUBROUTINE rhg_evp_rst 
    9301115 
     1116    
    9311117#else 
    9321118   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.