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 15531 – NEMO

Changeset 15531


Ignore:
Timestamp:
2021-11-23T19:09:14+01:00 (7 months ago)
Author:
vancop
Message:

Bugfixed VP rheology for trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_rhg_vp.F90

    r15519 r15531  
    4141   PUBLIC   ice_dyn_rhg_vp   ! called by icedyn_rhg.F90 
    4242 
    43  
    44    LOGICAL  ::   lp_zebra_vp =.TRUE.      ! activate zebra (solve the linear system problem every odd j-band, then one every even one) 
     43   INTEGER  ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp) 
     44   LOGICAL  ::   lp_zebra_vp =.TRUE. ! activate zebra (solve the linear system problem every odd j-band, then one every even one) 
    4545   REAL(wp) ::   zrelaxu_vp=0.95     ! U-relaxation factor (MV: can probably be merged with V-factor once ok) 
    4646   REAL(wp) ::   zrelaxv_vp=0.95     ! V-relaxation factor  
    4747   REAL(wp) ::   zuerr_max_vp=0.80   ! maximum velocity error, above which a forcing error is considered and solver is stopped 
    48    REAL(wp) ::   zuerr_min_vp=1.e-04   ! minimum velocity error, beyond which convergence is assumed 
     48   REAL(wp) ::   zuerr_min_vp=1.e-04 ! minimum velocity error, beyond which convergence is assumed 
    4949 
    5050   !! for convergence tests 
    5151   INTEGER ::   ncvgid        ! netcdf file id 
    52    INTEGER ::   nvarid_ures 
    53    INTEGER ::   nvarid_vres 
    54    INTEGER ::   nvarid_velres 
    55    INTEGER ::   nvarid_udif 
    56    INTEGER ::   nvarid_vdif 
    57    INTEGER ::   nvarid_veldif 
     52   INTEGER ::   nvarid_ures, nvarid_vres, nvarid_velres 
     53   INTEGER ::   nvarid_uerr_max, nvarid_verr_max, nvarid_velerr_max 
     54   INTEGER ::   nvarid_umad, nvarid_vmad, nvarid_velmad 
     55   INTEGER ::   nvarid_umad_outer, nvarid_vmad_outer, nvarid_velmad_outer 
    5856   INTEGER ::   nvarid_mke 
    59    INTEGER ::   nvarid_ures_xy, nvarid_vres_xy 
    6057 
    6158   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   fimask   ! mask at F points for the ice 
     
    8885      !! 
    8986      !!  f1(u) = g1(v) 
    90       !!  f2(v) = g2(v) 
     87      !!  f2(v) = g2(u) 
    9188      !! 
    9289      !!  The right-hand side (RHS) is explicit 
     
    141138      ! 
    142139      INTEGER ::   ji, ji2, jj, jj2, jn          ! dummy loop indices 
    143       INTEGER ::   jter, i_out, i_inn  !  
     140      INTEGER ::   i_out, i_inn, i_inn_tot  !  
    144141      INTEGER ::   ji_min, jj_min      ! 
    145142      INTEGER ::   nn_zebra_vp         ! number of zebra steps 
    146143 
    147       INTEGER ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp)       
    148144      ! 
    149145      REAL(wp) ::   zrhoco                                              ! rho0 * rn_cio 
     
    152148      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    153149      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass and volume 
    154       REAL(wp) ::   zdeltat, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars 
    155       REAL(wp) ::   zp_deltastar_f                                      !  
     150      REAL(wp) ::   zds2, zdt, zdt2, zdiv, zdiv2                        ! temporary scalars 
     151      REAL(wp) ::   zp_delstar_f                                        !  
    156152      REAL(wp) ::   zu_cV, zv_cU                                        !  
    157153      REAL(wp) ::   zfac, zfac1, zfac2, zfac3 
     
    164160      REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step 
    165161      ! 
    166       REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points 
    167       REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! Tension 
    168       REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta* at T points 
     162      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltat, zdelstar_t             ! Delta & Delta* at T-points 
     163      REAL(wp), DIMENSION(jpi,jpj) ::   ztens, zshear                   ! Tension, shear 
     164      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delstar_t                    ! P/delta* at T points 
    169165      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zet                        ! Viscosity pre-factors at T points 
    170166      REAL(wp), DIMENSION(jpi,jpj) ::   zef                             ! Viscosity pre-factor at F point 
     
    194190      REAL(wp), DIMENSION(jpi,jpj) ::   zFU, zFU_prime, zBU_prime       ! Rearranged linear system coefficients, U equation 
    195191      REAL(wp), DIMENSION(jpi,jpj) ::   zFV, zFV_prime, zBV_prime       ! Rearranged linear system coefficients, V equation 
    196       REAL(wp), DIMENSION(jpi,jpj) ::   zCU_prime, zCV_prime            ! Rearranged linear system coefficients, V equation 
    197192!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast) 
    198193!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast) 
    199194     ! 
    200       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00, zmsk15 
     195      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00 
    201196      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! mask for lots of ice (1), little ice (0) 
    202197      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence (1), no ice (0) 
     
    206201      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    207202      !! --- diags 
    208       REAL(wp) ::   zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 
     203      REAL(wp)                     ::   zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 
    209204      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f           ! stress tensor components 
    210205      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p 
     
    214209      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp, zdiag_yatrp         ! X/Y-component of area transport (m2/s, SIMIP) 
    215210 
    216  
    217       CALL ctl_stop( 'STOP', 'icedyn_rhg_vp: stop because vp rheology is an ongoing work and should not be used' ) 
     211      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zvel_res                         ! Residual of the linear system at last iteration 
     212      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zvel_diff                        ! Absolute velocity difference @last outer iteration 
     213                         
    218214       
    219215      !!---------------------------------------------------------------------------------------------------------------------- 
    220 !!$      ! DEBUG put all forcing terms to zero 
    221 !!$         ! air-ice drag 
    222 !!$         utau_ice(:,:) = 0._wp 
    223 !!$         vtau_ice(:,:) = 0._wp 
    224 !!$         ! coriolis 
    225 !!$         ff_t(:,:) = 0._wp 
    226 !!$         ! ice-ocean drag 
    227 !!$         rn_cio = 0._wp 
    228 !!$         ! ssh  
    229 !!$         ! done line 330 !!! dont forget to act there 
    230 !!$      ! END DEBUG 
    231216 
    232217      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' 
     
    243228         zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice 
    244229      END_2D 
    245       IF( nn_rhg_chkcvg > 0 ) THEN 
    246          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    247             zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 
    248          END_2D 
    249       ENDIF 
    250230       
    251231      IF ( lp_zebra_vp ) THEN; nn_zebra_vp = 2 
     
    267247      IF( nn_rhg_chkcvg /= 0 ) THEN 
    268248       
    269          ! ice area for global mean kinetic energy 
    270          zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) ) ! global ice area (km2) 
     249         ! ice area for global mean kinetic energy (m2) 
     250         zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 
    271251          
    272252      ENDIF 
     
    279259 
    280260      zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp; zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp 
    281       zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp; 
    282       zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp; 
    283       zrhsu(:,:) = 0._wp; zrhsv(:,:) = 0._wp 
    284       zf_rhsu(:,:) = 0._wp; zf_rhsv(:,:) = 0._wp 
     261      zrhsu  (:,:)  = 0._wp; zrhsv  (:,:)  = 0._wp; zf_rhsu(:,:)  = 0._wp; zf_rhsv(:,:)  = 0._wp 
     262      zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp 
     263      zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp 
    285264 
    286265      !------------------------------------------------------------------------------! 
     
    292271      CALL ice_strength ! strength at T points 
    293272       
    294       !------------------------------ 
    295       ! -- F-mask       (code from EVP) 
    296       !------------------------------ 
     273      !--------------------------- 
     274      ! -- F-mask (code from EVP) 
     275      !--------------------------- 
    297276      IF( kt == nit000 ) THEN 
    298277         ! MartinV:  
     
    328307      !    embedded sea ice: compute representative ice top surface 
    329308      !    non-embedded sea ice: use ocean surface for slope calculation 
    330       zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
    331 !!$      zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!! 
     309      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)     
    332310 
    333311      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    334          zm1          = ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )  ! Ice/snow mass at U-V points 
    335          zmf  (ji,jj) = zm1 * ff_t(ji,jj)                            ! Coriolis at T points (m*f) 
     312         zmt(ji,jj) = rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj)       ! Snow and ice mass at T-point 
     313         zmf(ji,jj) = zmt(ji,jj) * ff_t(ji,jj)                      ! Coriolis factor at T points (m*f) 
    336314      END_2D 
    337315       
    338       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     316      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
    339317 
    340318         ! Ice fraction at U-V points 
    341319         za_iU(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    342320         za_iV(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    343           
    344          ! Ice/snow mass at U-V points 
    345          zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) ) 
    346          zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) ) 
    347          zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )          
     321 
     322         ! Snow and ice mass at U-V points 
     323         zm1             = zmt(ji,jj) 
     324         zm2             = zmt(ji+1,jj) 
     325         zm3             = zmt(ji,jj+1) 
    348326         zmassU          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    349327         zmassV          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     
    379357         ELSE                                                      ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
    380358 
    381 ! MV TEST DEBUG 
    382 !!$            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji+1,jj)  <= zmmin )     .AND.  & 
    383 !!$               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji+1,jj) <= zamin ) )    THEN   ;   zmsk01x(ji,jj) = 0._wp 
    384 !!$            ELSE                                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
    385 !!$ 
    386 !!$            IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji,jj+1)  <= zmmin )     .AND.  & 
    387 !!$               & ( at_i(ji,jj)  <= zamin .OR. at_i(ji,jj+1) <= zamin ) )    THEN   ;   zmsk01y(ji,jj) = 0._wp 
    388 !!$            ELSE                                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
    389 ! END MV TEST DEBUG 
    390  
    391       END_2D 
    392  
    393 !!$      CALL iom_put( 'zmsk00x'    , zmsk00x  )   ! MV DEBUG 
    394 !!$      CALL iom_put( 'zmsk00y'    , zmsk00y  )   ! MV DEBUG 
    395 !!$      CALL iom_put( 'zmsk01x'    , zmsk01x  )   ! MV DEBUG 
    396 !!$      CALL iom_put( 'zmsk01y'    , zmsk01y  )   ! MV DEBUG 
    397 !!$      CALL iom_put( 'ztaux_ai'   , ztaux_ai )   ! MV DEBUG 
    398 !!$      CALL iom_put( 'ztauy_ai'   , ztauy_ai )   ! MV DEBUG 
    399 !!$      CALL iom_put( 'zspgU'      , zspgU    )   ! MV DEBUG 
    400 !!$      CALL iom_put( 'zspgV'      , zspgV    )   ! MV DEBUG 
    401 !!$             
     359      END_2D   
     360             
    402361      !------------------------------------------------------------------------------! 
    403362      ! 
     
    409368      zv_c(:,:) = v_ice(:,:) 
    410369       
    411       jter = 0 
     370      i_inn_tot = 0 
    412371 
    413372      DO i_out = 1, nn_vp_nout 
    414373 
    415          IF( lwp )   WRITE(numout,*) ' outer loop  i_out : ', i_out 
    416                 
    417374         ! Velocities used in the non linear terms are the average of the past two iterates 
    418          ! u_it = 0.5 * ( u_{it-1} + u_{it-2}) 
     375         ! u_it = 0.5 * ( u_{it-1} + u_{it-2} ) 
    419376         ! Also used in Hibler and Ackley (1983); Zhang and Hibler (1997); Lemieux and Tremblay (2009) 
    420377         zu_c(:,:) = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) ) 
     
    428385         ! In the outer loop, one needs to update all RHS terms 
    429386         ! with explicit velocity dependencies (viscosities, coriolis, ocean stress) 
    430          ! as a function of uc 
    431          ! 
     387         ! as a function of "current" velocities (uc, vc) 
    432388       
    433389         !------------------------------------------ 
     
    436392 
    437393         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
    438          DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
    439  
     394         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! 1->jpi-1 
     395          
     396            ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
    440397            ! shear at F points 
    441398            zds(ji,jj) = ( ( zu_c(ji,jj+1) * r1_e1u(ji,jj+1) - zu_c(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     
    445402         END_2D 
    446403 
    447          CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan 
    448 !!$         CALL iom_put( 'zds'          , zds      )   ! MV DEBUG 
    449  
    450          IF( lwp )   WRITE(numout,*) ' outer loop  1a i_out : ', i_out 
    451  
    452          !DO jj = 2, jpj - 1    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
    453          !   DO ji = 2, jpi - 1 !  
    454  
    455 ! MV DEBUG 
    456          DO jj = 2, jpj        ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
    457             DO ji = 2, jpi     !  
    458 ! END MV DEBUG 
    459  
    460                ! shear**2 at T points (doc eq. A16) 
    461                zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
    462                   &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
    463                   &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
     404         CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! necessary, zds2 uses jpi/jpj values for zds  
     405 
     406         DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! 2 -> jpj; 2,jpi !!! CHECK !!! 
     407            ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
     408 
     409            ! shear**2 at T points (doc eq. A16) 
     410            zds2  = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     411               &    + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
     412               &    ) * 0.25_wp * r1_e1e2t(ji,jj) 
    464413               
    465                ! divergence at T points 
    466                zdiv  = ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj)   & 
    467                   &    + e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1)   & 
    468                   &    ) * r1_e1e2t(ji,jj) 
    469                zdiv2 = zdiv * zdiv 
     414            ! divergence at T points 
     415            zdiv  = ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj)   & 
     416               &    + e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1)   & 
     417               &    ) * r1_e1e2t(ji,jj) 
     418            zdiv2 = zdiv * zdiv 
    470419                
    471                ! tension at T points 
    472                zdt  = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    473                   &   - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    474                   &   ) * r1_e1e2t(ji,jj) 
    475                zdt2 = zdt * zdt 
     420            ! tension at T points 
     421            zdt   = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     422               &    - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     423               &    ) * r1_e1e2t(ji,jj) 
     424            zdt2 = zdt * zdt 
    476425                
    477                ! delta at T points 
    478                zdeltat = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     426            ! delta at T points 
     427            zdeltat(ji,jj)        = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
    479428                
    480                ! delta* at T points (following Lemieux and Dupont, GMD 2020) 
    481                zdeltastar_t(ji,jj) = zdeltat + rn_creepl 
     429            ! delta* at T points (following Lemieux and Dupont, GMD 2020) 
     430            zdelstar_t(ji,jj)     = zdeltat(ji,jj) + rn_creepl ! OPT zdelstar_t can be totally removed and put into next line directly. Could change results 
    482431               
    483                ! P/delta at T-points 
    484                zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltastar_t(ji,jj) 
     432            ! P/delta* at T-points 
     433            zp_delstar_t(ji,jj)   = strength(ji,jj) / zdelstar_t(ji,jj) 
    485434                
    486                ! Temporary zzt and zet factors at T-points 
    487                zzt(ji,jj)     = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj) 
    488                zet(ji,jj)     = zzt(ji,jj)     * z1_ecc2  
     435            ! Temporary zzt and zet factors at T-points 
     436            zzt(ji,jj)            = zp_delstar_t(ji,jj) * r1_e1e2t(ji,jj) 
     437            zet(ji,jj)            = zzt(ji,jj)     * z1_ecc2  
    489438                           
    490             END DO 
    491          END DO 
    492           
    493          CALL lbc_lnk( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 
    494  
    495 !!$         CALL iom_put( 'zzt'        , zzt      )   ! MV DEBUG 
    496 !!$         CALL iom_put( 'zet'        , zet      )   ! MV DEBUG 
    497 !!$         CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG 
    498  
    499          IF( lwp )   WRITE(numout,*) ' outer loop  1b i_out : ', i_out 
    500  
    501          DO jj = 1, jpj - 1 
    502             DO ji = 1, jpi - 1 
    503                 
     439         END_2D 
     440          
     441         CALL lbc_lnk( 'icedyn_rhg_vp', zp_delstar_t , 'T', 1. ) ! necessary, used for ji = 1 and jj = 1 
     442 
     443         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )! 1-> jpj-1; 1->jpi-1 
     444          
    504445               ! P/delta* at F points 
    505                zp_deltastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) ) 
     446               zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) ) 
    506447                
    507448               ! Temporary zef factor at F-point 
    508                zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * fimask(ji,jj) 
    509  
    510             END DO 
    511          END DO 
    512           
    513          CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 
    514 !!$         CALL iom_put( 'zef'          , zef            ) ! MV DEBUG 
    515          IF( lwp )   WRITE(numout,*) ' outer loop  1c i_out : ', i_out 
    516  
     449               zef(ji,jj)      = zp_delstar_f * r1_e1e2f(ji,jj) * z1_ecc2 * fimask(ji,jj) * 0.5_wp 
     450                
     451         END_2D 
     452          
    517453         !--------------------------------------------------- 
    518454         ! -- Ocean-ice drag and Coriolis RHS contributions 
    519455         !--------------------------------------------------- 
    520456 
    521          DO jj = 2, jpj - 1 
    522              DO ji = 2, jpi - 1 
     457         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     458          
     459            !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 
     460            zu_cV            = 0.25_wp * ( zu_c(ji,jj) + zu_c(ji-1,jj) + zu_c(ji,jj+1) + zu_c(ji-1,jj+1) ) * vmask(ji,jj,1) 
     461            zv_cU            = 0.25_wp * ( zv_c(ji,jj) + zv_c(ji,jj-1) + zv_c(ji+1,jj) + zv_c(ji+1,jj-1) ) * umask(ji,jj,1) 
    523462                 
    524                 !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 
    525                 zu_cV            = 0.25_wp * ( zu_c(ji,jj) + zu_c(ji-1,jj) + zu_c(ji,jj+1) + zu_c(ji-1,jj+1) ) * vmask(ji,jj,1) 
    526                 zv_cU            = 0.25_wp * ( zv_c(ji,jj) + zv_c(ji,jj-1) + zv_c(ji+1,jj) + zv_c(ji+1,jj-1) ) * umask(ji,jj,1) 
     463            !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 
     464            zCwU(ji,jj)          = za_iU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) )  & 
     465              &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 
     466            zCwV(ji,jj)          = za_iV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) )  & 
     467              &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 
     468                  
     469            !--- Ocean-ice drag contributions to RHS  
     470            ztaux_oi_rhsu(ji,jj) = zCwU(ji,jj) * u_oce(ji,jj) 
     471            ztauy_oi_rhsv(ji,jj) = zCwV(ji,jj) * v_oce(ji,jj) 
    527472                 
    528                 !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 
    529                 zCwU(ji,jj)          = za_iU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) )  & 
    530                   &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 
    531                 zCwV(ji,jj)          = za_iV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) )  & 
    532                   &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 
    533                   
    534                 !--- Ocean-ice drag contributions to RHS  
    535                 ztaux_oi_rhsu(ji,jj) = zCwU(ji,jj) * u_oce(ji,jj) 
    536                 ztauy_oi_rhsv(ji,jj) = zCwV(ji,jj) * v_oce(ji,jj) 
    537                  
    538                 ! --- U-component of Coriolis Force (energy conserving formulation) 
    539                 ! Note Lemieux et al 2008 recommend to do that implicitly, but I don't really see how this could be done 
    540                 zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  & 
    541                            &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * zv_c(ji  ,jj) + e1v(ji  ,jj-1) * zv_c(ji  ,jj-1) )  & 
    542                            &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_c(ji+1,jj) + e1v(ji+1,jj-1) * zv_c(ji+1,jj-1) ) ) 
     473            !--- U-component of Coriolis Force (energy conserving formulation) 
     474            zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  & 
     475                       &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * zv_c(ji  ,jj) + e1v(ji  ,jj-1) * zv_c(ji  ,jj-1) )  & 
     476                       &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_c(ji+1,jj) + e1v(ji+1,jj-1) * zv_c(ji+1,jj-1) ) ) 
    543477                            
    544                 ! --- V-component of Coriolis Force (energy conserving formulation) 
    545                 zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  & 
    546                            &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * zu_c(ji,jj  ) + e2u(ji-1,jj  ) * zu_c(ji-1,jj  ) )  & 
    547                            &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) ) 
    548           
    549              END DO 
    550          END DO 
    551  
    552          IF( lwp )   WRITE(numout,*) ' outer loop  1d i_out : ', i_out 
    553           
    554          CALL lbc_lnk( 'icedyn_rhg_vp', zCwU ,  'U', -1., zCwV, 'V', -1. ) 
    555          CALL lbc_lnk( 'icedyn_rhg_vp', zCorU,  'U', -1., zCorV, 'V', -1. ) 
    556  
    557 !!$         CALL iom_put( 'zCwU'          , zCwU           ) ! MV DEBUG 
    558 !!$         CALL iom_put( 'zCwV'          , zCwV           ) ! MV DEBUG 
    559 !!$         CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG 
    560 !!$         CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG 
    561  
    562          IF( lwp )   WRITE(numout,*) ' outer loop  1f i_out : ', i_out 
    563           
    564          ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system,  
    565          ! so there is no need for lbclnk on drag and coriolis 
    566  
     478            !--- V-component of Coriolis Force (energy conserving formulation) 
     479            zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  & 
     480                       &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * zu_c(ji,jj  ) + e2u(ji-1,jj  ) * zu_c(ji-1,jj  ) )  & 
     481                       &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) ) 
     482 
     483         END_2D 
     484          
    567485         !------------------------------------- 
    568486         ! -- Internal stress RHS contribution 
    569487         !------------------------------------- 
    570488 
    571          ! --- Stress contributions at T-points          
    572          DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
    573             DO ji = 2, jpi !  
     489         ! --- Stress contributions at T-points 
     490         DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! 2 -> jpj; 2,jpi !!! CHECK !!! 
     491          
     492         ! loop to jpi,jpj to avoid making a communication for zs1 & zs2 
    574493             
    575                ! sig1 contribution to RHS of U-equation at T-points 
    576                zs1_rhsu(ji,jj) =   zzt(ji,jj) * ( e1v(ji,jj)    * zv_c(ji,jj) - e1v(ji,jj-1)    * zv_c(ji,jj-1) - 1.0_wp ) 
     494            ! sig1 contribution to RHS of U-equation at T-points 
     495            zs1_rhsu(ji,jj) =   zzt(ji,jj) * ( e1v(ji,jj)    * zv_c(ji,jj) - e1v(ji,jj-1)    * zv_c(ji,jj-1) )   & 
     496                            &                - zp_delstar_t(ji,jj) * zdeltat(ji,jj) 
    577497                                             
    578                ! sig2 contribution to RHS of U-equation at T-points             
    579                zs2_rhsu(ji,jj) = - zet(ji,jj) * ( r1_e1v(ji,jj) * zv_c(ji,jj) - r1_e1v(ji,jj-1) * zv_c(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  
    580  
    581                ! sig1 contribution to RHS of V-equation at T-points 
    582                zs1_rhsv(ji,jj) =   zzt(ji,jj) * ( e2u(ji,jj)    * zu_c(ji,jj) - e2u(ji-1,jj)    * zu_c(ji-1,jj) - 1.0_wp ) 
    583  
    584                ! sig2 contribution to RHS of V-equation  at T-points 
    585                zs2_rhsv(ji,jj) =   zet(ji,jj) * ( r1_e2u(ji,jj) * zu_c(ji,jj) - r1_e2u(ji-1,jj) * zu_c(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) 
    586  
    587             END DO 
    588          END DO 
    589  
    590 !!$         CALL iom_put( 'zs1_rhsu'      , zs1_rhsu       ) ! MV DEBUG 
    591 !!$         CALL iom_put( 'zs2_rhsu'      , zs2_rhsu       ) ! MV DEBUG 
    592 !!$         CALL iom_put( 'zs1_rhsv'      , zs1_rhsv       ) ! MV DEBUG 
    593 !!$         CALL iom_put( 'zs2_rhsv'      , zs2_rhsv       ) ! MV DEBUG 
    594           
    595          ! a priori, no lbclnk, because rhsu is only used in the inner domain 
    596           
    597          ! --- Stress contributions at f-points          
     498            ! sig2 contribution to RHS of U-equation at T-points             
     499            zs2_rhsu(ji,jj) = - zet(ji,jj) * ( r1_e1v(ji,jj) * zv_c(ji,jj) - r1_e1v(ji,jj-1) * zv_c(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  
     500 
     501            ! sig1 contribution to RHS of V-equation at T-points 
     502            zs1_rhsv(ji,jj) =   zzt(ji,jj) * ( e2u(ji,jj)    * zu_c(ji,jj) - e2u(ji-1,jj)    * zu_c(ji-1,jj) )   &  
     503                            &                - zp_delstar_t(ji,jj) * zdeltat(ji,jj) 
     504 
     505            ! sig2 contribution to RHS of V-equation  at T-points 
     506            zs2_rhsv(ji,jj) =   zet(ji,jj) * ( r1_e2u(ji,jj) * zu_c(ji,jj) - r1_e2u(ji-1,jj) * zu_c(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) 
     507 
     508         END_2D 
     509                   
     510         ! --- Stress contributions at F-points          
    598511         ! MV NOTE: I applied fimask on zds, by mimetism on EVP, but without deep understanding of what I was doing 
    599512         ! My guess is that this is the way to enforce boundary conditions on strain rate tensor 
    600513 
    601          IF( lwp )   WRITE(numout,*) ' outer loop 2 i_out : ', i_out 
    602  
    603          DO jj = 1, jpj - 1 
    604             DO ji = 1, jpi - 1 
    605                 
    606                ! sig12 contribution to RHS of U equation at F-points  
    607                zs12_rhsu(ji,jj) = - zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) - r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * fimask(ji,jj) 
    608                 
    609                ! sig12 contribution to RHS of V equation at F-points 
    610                zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) - r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * fimask(ji,jj) 
    611  
    612             END DO 
    613          END DO 
    614  
    615          CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) 
    616          CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. ) 
    617  
    618 !!$         CALL iom_put( 'zs12_rhsu'     , zs12_rhsu      ) ! MV DEBUG 
    619 !!$         CALL iom_put( 'zs12_rhsv'     , zs12_rhsv      ) ! MV DEBUG 
    620  
    621          ! a priori, no lbclnk, because rhsu are only used in the inner domain 
    622  
     514         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! 1->jpi-1 
     515          
     516            ! sig12 contribution to RHS of U equation at F-points  
     517            zs12_rhsu(ji,jj) =   zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) + r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * fimask(ji,jj) 
     518 
     519            ! sig12 contribution to RHS of V equation at F-points 
     520            zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) + r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * fimask(ji,jj) 
     521 
     522         END_2D 
     523          
    623524         ! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002) 
    624525         ! OPT: merge with next loop and use intermediate scalars for zf_rhsu 
    625           
    626          DO jj = 2, jpj - 1 
    627             DO ji = 2, jpi - 1                
     526         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     527          
    628528               ! --- U component of internal force contribution to RHS at U points 
    629529               zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * &  
     
    635535               zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * & 
    636536                  &           (    e1v(ji,jj)    * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) )                                                                 & 
    637                   &      +         r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj) * e1t(ji,jj)     * zs2_rhsv(ji,jj) )     & 
     537                  &      -         r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj)   * e1t(ji,jj)   * zs2_rhsv(ji,jj) )     & 
    638538                  &      + 2._wp * r1_e2v(ji,jj) * ( e2f(ji,jj)   * e2f(ji,jj)   * zs12_rhsv(ji,jj)  - e2f(ji-1,jj) * e2f(ji-1,jj) * zs12_rhsv(ji-1,jj) ) ) 
    639                    
    640             END DO 
    641          END DO 
    642  
    643 !!$         CALL iom_put( 'zf_rhsu'       , zf_rhsu        ) ! MV DEBUG 
    644 !!$         CALL iom_put( 'zf_rhsv'       , zf_rhsv        ) ! MV DEBUG 
     539 
     540         END_2D 
    645541          
    646542         !--------------------------- 
     
    649545         ! 
    650546         ! OPT: could use intermediate scalars to reduce memory access 
    651          DO jj = 2, jpj - 1 
    652             DO ji = 2, jpi - 1 
    653              
    654                ! still miss ice ocean stress and acceleration contribution 
    655                zrhsu(ji,jj) = zmU_t(ji,jj) + ztaux_ai(ji,jj) + ztaux_oi_rhsu(ji,jj) + zspgU(ji,jj) + zCorU(ji,jj) + zf_rhsu(ji,jj) 
    656                zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsu(ji,jj) 
    657  
    658             END DO 
    659          END DO 
    660           
    661          CALL lbc_lnk( 'icedyn_rhg_vp', zrhsu, 'U', -1., zrhsv, 'V',  -1.) 
    662          CALL lbc_lnk( 'icedyn_rhg_vp', zmU_t, 'U', -1., zmV_t, 'V',  -1.) 
    663          CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V',  -1.) 
    664  
    665 !!$         CALL iom_put( 'zmU_t'         , zmU_t          ) ! MV DEBUG 
    666 !!$         CALL iom_put( 'zmV_t'         , zmV_t          ) ! MV DEBUG 
    667 !!$         CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu  ) ! MV DEBUG 
    668 !!$         CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv  ) ! MV DEBUG 
    669 !!$         CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG 
    670 !!$         CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG 
    671 !!$          
    672          ! inner domain calculations -> no lbclnk 
    673  
    674          IF( lwp )   WRITE(numout,*) ' outer loop 4 i_out : ', i_out 
    675       
     547         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     548                      
     549            zrhsu(ji,jj) = zmU_t(ji,jj) + ztaux_ai(ji,jj) + ztaux_oi_rhsu(ji,jj) + zspgU(ji,jj) + zCorU(ji,jj) + zf_rhsu(ji,jj) 
     550            zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsv(ji,jj) 
     551 
     552         END_2D 
     553          
    676554         !---------------------------------------------------------------------------------------! 
    677555         ! 
     
    691569         !         only zzt and zet are iteration-dependent, other only depend on scale factors 
    692570                   
    693          DO ji = 2, jpi - 1 ! internal domain do loop 
    694             DO jj = 2, jpj - 1 
    695  
    696                !------------------------------------- 
    697                ! -- Internal forces LHS contribution 
    698                !------------------------------------- 
    699                ! 
    700                ! --- U-component 
    701                ! 
    702                ! "T" factors (intermediate results) 
    703                ! 
    704                zfac       = 0.5_wp * r1_e1e2u(ji,jj) 
    705                zfac1      =         zfac * e2u(ji,jj) 
    706                zfac2      =         zfac * r1_e2u(ji,jj) 
    707                zfac3      = 2._wp * zfac * r1_e1u(ji,jj) 
    708  
    709                zt12U      = - zfac1 * zzt(ji+1,jj) 
    710                zt11U      =   zfac1 * zzt(ji,jj) 
    711           
    712                zt22U      = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 
    713                zt21U      =   zfac2 * zet(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj) 
    714           
    715                zt122U     = - zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj) 
    716                zt121U     =   zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 
     571         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     572          
     573            !------------------------------------- 
     574            ! -- Internal forces LHS contribution 
     575            !------------------------------------- 
     576            ! 
     577            ! --- U-component 
     578            ! 
     579            ! "T" factors (intermediate results) 
     580            ! 
     581            zfac       = 0.5_wp * r1_e1e2u(ji,jj) 
     582            zfac1      =         zfac * e2u(ji,jj) 
     583            zfac2      =         zfac * r1_e2u(ji,jj) 
     584            zfac3      = 2._wp * zfac * r1_e1u(ji,jj) 
     585 
     586            zt11U      =   zfac1 * zzt(ji,jj) 
     587            zt12U      =   zfac1 * zzt(ji+1,jj) 
     588          
     589            zt21U      =   zfac2 * zet(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj) 
     590            zt22U      =   zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 
     591          
     592            zt121U     =   zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 
     593            zt122U     =   zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj) 
    717594                
    718                ! 
    719                ! Linear system coefficients 
    720                ! 
    721                zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj) 
    722                zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj) 
    723                zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj) 
    724           
    725                zDU(ji,jj) =   zt121U * r1_e1u(ji,jj-1) 
    726                zEU(ji,jj) =   zt122U * r1_e1u(ji,jj+1) 
     595            ! 
     596            ! Linear system coefficients 
     597            ! 
     598            zAU(ji,jj) = -   zt11U           * e2u(ji-1,jj) -   zt21U          * r1_e2u(ji-1,jj) 
     599            zBU(ji,jj) =   ( zt11U + zt12U ) * e2u(ji,jj)   + ( zt21U + zt22U ) * r1_e2u(ji,jj)   + ( zt121U + zt122U ) * r1_e1u(ji,jj) 
     600            zCU(ji,jj) = -   zt12U           * e2u(ji+1,jj) -   zt22U          * r1_e2u(ji+1,jj) 
     601          
     602            zDU(ji,jj) =     zt121U * r1_e1u(ji,jj-1) 
     603            zEU(ji,jj) =     zt122U * r1_e1u(ji,jj+1) 
    727604               
    728                ! 
    729                ! --- V-component 
    730                ! 
    731                ! "T" factors (intermediate results) 
    732                ! 
    733                zfac       = 0.5_wp * r1_e1e2v(ji,jj) 
    734                zfac1      =         zfac * e2v(ji,jj) 
    735                zfac2      =         zfac * r1_e1v(ji,jj) 
    736                zfac3      = 2._wp * zfac * r1_e2v(ji,jj) 
    737           
    738                zt12V      = - zfac1 * zzt(ji,jj+1) 
    739                zt11V      =   zfac1 * zzt(ji,jj) 
    740           
    741                zt22V      =   zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) 
    742                zt21V      = - zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj) 
    743           
    744                zt122V     =   zfac3 * zef(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj) 
    745                zt121V     = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 
    746           
    747                ! 
    748                ! Linear system coefficients 
    749                ! 
    750                zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1) 
    751                zBV(ji,jj) =  ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj) 
    752                zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1) 
    753           
    754                zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here 
    755                zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj) 
     605            ! 
     606            ! --- V-component 
     607            ! 
     608            ! "T" factors (intermediate results) 
     609            ! 
     610            zfac       = 0.5_wp * r1_e1e2v(ji,jj) 
     611            zfac1      =         zfac * e1v(ji,jj) 
     612            zfac2      =         zfac * r1_e1v(ji,jj) 
     613            zfac3      = 2._wp * zfac * r1_e2v(ji,jj) 
     614 
     615            zt11V      =   zfac1 * zzt(ji,jj) 
     616            zt12V      =   zfac1 * zzt(ji,jj+1) 
     617 
     618            zt21V      =   zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj) 
     619            zt22V      =   zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) 
     620          
     621            zt121V     =   zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 
     622            zt122V     =   zfac3 * zef(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj) 
     623 
     624            ! 
     625            ! Linear system coefficients 
     626            ! 
     627            zAV(ji,jj) = -   zt11V           * e1v(ji,jj-1) -   zt21V          * r1_e1v(ji,jj-1) 
     628            zBV(ji,jj) =   ( zt11V + zt12V ) * e1v(ji,jj)   + ( zt21V + zt22V ) * r1_e1v(ji,jj)   + ( zt122V + zt121V ) * r1_e2v(ji,jj) 
     629            zCV(ji,jj) = -   zt12V           * e1v(ji,jj+1) -   zt22V          * r1_e1v(ji,jj+1) 
     630 
     631            zDV(ji,jj) =     zt121V * r1_e2v(ji-1,jj) 
     632            zEV(ji,jj) =    zt122V * r1_e2v(ji+1,jj) 
    756633                   
    757                !----------------------------------------------------- 
    758                ! -- Ocean-ice drag and acceleration LHS contribution 
    759                !----------------------------------------------------- 
    760                zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj) 
    761                zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 
    762           
    763             END DO 
    764          END DO 
    765  
    766          CALL lbc_lnk( 'icedyn_rhg_vp', zAU  , 'U', 1., zAV  , 'V',  1. ) 
    767          CALL lbc_lnk( 'icedyn_rhg_vp', zBU  , 'U', 1., zBV  , 'V',  1. ) 
    768          CALL lbc_lnk( 'icedyn_rhg_vp', zCU  , 'U', 1., zCV  , 'V',  1. ) 
    769          CALL lbc_lnk( 'icedyn_rhg_vp', zDU  , 'U', 1., zDV  , 'V',  1. ) 
    770          CALL lbc_lnk( 'icedyn_rhg_vp', zEU  , 'U', 1., zEV  , 'V',  1. ) 
    771                 
    772 !!$         CALL iom_put( 'zAU'           , zAU            ) ! MV DEBUG 
    773 !!$         CALL iom_put( 'zBU'           , zBU            ) ! MV DEBUG 
    774 !!$         CALL iom_put( 'zCU'           , zCU            ) ! MV DEBUG 
    775 !!$         CALL iom_put( 'zDU'           , zDU            ) ! MV DEBUG 
    776 !!$         CALL iom_put( 'zEU'           , zEU            ) ! MV DEBUG 
    777 !!$         CALL iom_put( 'zAV'           , zAV            ) ! MV DEBUG 
    778 !!$         CALL iom_put( 'zBV'           , zBV            ) ! MV DEBUG 
    779 !!$         CALL iom_put( 'zCV'           , zCV            ) ! MV DEBUG 
    780 !!$         CALL iom_put( 'zDV'           , zDV            ) ! MV DEBUG 
    781 !!$         CALL iom_put( 'zEV'           , zEV            ) ! MV DEBUG 
    782  
     634            !----------------------------------------------------- 
     635            ! -- Ocean-ice drag and acceleration LHS contribution 
     636            !----------------------------------------------------- 
     637            zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj) 
     638            zBV(ji,jj) = zBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 
     639          
     640         END_2D 
     641          
    783642      !------------------------------------------------------------------------------! 
    784643      ! 
     
    793652         DO i_inn = 1, nn_vp_ninn ! inner loop iterations 
    794653 
    795             IF( lwp )   WRITE(numout,*) ' inner loop 1 i_inn : ', i_inn 
    796           
    797654            !--- mitgcm computes initial value of residual here... 
    798655 
    799             jter             = jter + 1 
    800             ! l_full_nf_update = jter == nn_nvp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
    801  
    802             zu_b(:,:)      = u_ice(:,:) ! velocity at previous sub-iterate 
    803             zv_b(:,:)      = v_ice(:,:) 
    804  
    805 !           zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp 
    806 !           zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp 
     656            i_inn_tot             = i_inn_tot + 1 
     657            ! l_full_nf_update = i_inn_tot == nn_nvp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
     658 
     659            zu_b(:,:)       = u_ice(:,:) ! velocity at previous inner-iterate 
     660            zv_b(:,:)       = v_ice(:,:) 
    807661 
    808662            IF ( ll_u_iterate .OR. ll_v_iterate )   THEN 
     
    817671                  ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F 
    818672 
    819                   zFU(:,:)       = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; zCU_prime(:,:) = 0._wp 
     673                  zFU(:,:)       = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp;  
    820674                   
    821675                  DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) 
     
    826680                     ELSE                  ;   jj_min = 3 
    827681                     ENDIF 
    828  
    829                      IF ( lwp ) WRITE(numout,*) ' Into the U-zebra loop at step jn = ', jn, ', with jj_min = ', jj_min 
    830682 
    831683                     DO jj = jj_min, jpj - 1, nn_zebra_vp 
     
    835687                        !------------------------ 
    836688                        DO ji = 2, jpi - 1     
    837  
    838                            ! boundary condition substitution 
     689                           ! note: these are key lines linking information between processors 
     690                           ! u_ice/v_ice need to be lbc_linked 
     691 
     692                           ! sub-domain boundary condition substitution 
    839693                           ! see Zhang and Hibler, 1997, Appendix B 
    840694                           zAA3 = 0._wp 
     
    852706                     END DO 
    853707                      
    854                      CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U',  1. ) 
    855                       
    856708                     !--------------- 
    857709                     ! Forward sweep 
     
    859711                     DO jj = jj_min, jpj - 1, nn_zebra_vp 
    860712       
     713                        zBU_prime(2,jj)     = zBU(2,jj) 
     714                        zFU_prime(2,jj)     = zFU(2,jj) 
     715 
    861716                        DO ji = 3, jpi - 1 
    862717 
     
    869724 
    870725                     END DO 
    871  
    872                      CALL lbc_lnk( 'icedyn_rhg_vp', zFU_prime, 'U',  1., zBU_prime, 'U', 1. ) 
    873   
     726                                                                                                      
    874727                     !----------------------------- 
    875728                     ! Backward sweep & relaxation 
     
    879732                     
    880733                        ! --- Backward sweep  
     734 
    881735                        ! last row  
    882736                        zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) ) 
    883737                        u_ice(jpi-1,jj)    = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) &  
    884738                                           &            * umask(jpi-1,jj,1) 
    885                         DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop 
     739 
     740                        DO ji = jpi - 2 , 2, -1 ! all other rows    !  ---> original backward loop 
    886741                           zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 
    887742                           u_ice(ji,jj)    = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1)   &  
     
    889744                        END DO 
    890745 
    891                         !--- Relaxation 
    892                         ! and velocity masking for little-ice and no-ice cases 
     746                        !--- Relaxation and masking (for low-ice/no-ice cases) 
    893747                        DO ji = 2, jpi - 1     
    894748                         
     
    902756 
    903757                     END DO ! jj 
     758 
     759                     CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1. ) 
    904760                      
    905761                  END DO ! zebra loop 
     
    917773                  !!! ZH97 explain it is critical for convergence speed 
    918774 
    919                   zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; zCV_prime(:,:) = 0._wp 
     775                  zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp;  
    920776 
    921777                  DO jn = 1, nn_zebra_vp ! "zebra" loop 
     
    925781                     ENDIF 
    926782 
    927                      IF ( lwp ) WRITE(numout,*) ' Into the V-zebra loop at step jn = ', jn, ', with ji_min = ', ji_min 
    928           
    929783                     DO ji = ji_min, jpi - 1, nn_zebra_vp  
    930784                      
     
    934788                        DO jj = 2, jpj - 1 
    935789 
    936                            ! boundary condition substitution (check it is correctly applied !!!) 
     790                           ! subdomain boundary condition substitution (check it is correctly applied !!!) 
    937791                           ! see Zhang and Hibler, 1997, Appendix B 
    938792                           zAA3 = 0._wp 
     
    950804                     END DO 
    951805 
    952                      CALL lbc_lnk( 'icedyn_rhg_vp', zFV, 'V',  1.) 
    953                       
    954806                     !--------------- 
    955807                     ! Forward sweep 
     
    957809                     DO ji = ji_min, jpi - 1, nn_zebra_vp  
    958810                      
    959                         DO jj = 3, jpj - 1 
     811                        zBV_prime(ji,2)     = zBV(ji,2) 
     812                        zFV_prime(ji,2)     = zFV(ji,2) 
     813 
     814                        DO jj = 3, jpj - 1  
    960815 
    961816                           zfac             = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) 
     
    968823                     END DO 
    969824 
    970                      CALL lbc_lnk( 'icedyn_rhg_vp', zFV_prime, 'V',  1., zBV_prime, 'V', 1. ) 
    971                       
    972825                     !----------------------------- 
    973826                     ! Backward sweep & relaxation 
     
    988841                        END DO             
    989842                                                    
    990                         ! --- Relaxation  & masking (should it be now or later) 
     843                        ! --- Relaxation & masking  
    991844                        DO jj = 2, jpj - 1 
    992845                         
     
    1000853                         
    1001854                     END DO ! ji 
     855 
     856                     CALL lbc_lnk( 'icedyn_rhg_vp', v_ice, 'V', -1. ) 
    1002857                      
    1003858                  END DO ! zebra loop 
     
    1005860               ENDIF !   ll_v_iterate 
    1006861 
    1007                CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
     862               ! I suspect the communication should go into the zebra loop if we want reproducibility 
    1008863                               
    1009864               !-------------------------------------------------------------------------------------- 
     
    1016871               ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check 
    1017872               ! then we must optimize the use of the mpp_max, which is prohibitive                             
    1018                zuerr_max = 0._wp 
     873               zuerr_max  = 0._wp 
    1019874                                
    1020875               IF ( ll_u_iterate .AND. MOD ( i_inn, nn_vp_chkcvg ) == 0 ) THEN 
     
    1022877                  ! - Maximum U-velocity difference                
    1023878                  zuerr(:,:) = 0._wp 
    1024                   DO jj = 2, jpj - 1 
    1025                      DO ji = 2, jpi - 1 
    1026                         zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 
    1027                      END DO 
    1028                   END DO 
     879                  DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     880                   
     881                     zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1)  
     882                   
     883                  END_2D 
     884          
    1029885                  zuerr_max = MAXVAL( zuerr ) 
    1030886                  CALL mpp_max( 'icedyn_rhg_evp', zuerr_max )   ! max over the global domain - damned! 
    1031                    
    1032                   ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 
     887 
     888                  ! - Stop if max error is too large ("safeguard against bad forcing" of original Zhang routine) 
    1033889                  IF ( i_inn > 1 .AND. zuerr_max > zuerr_max_vp ) THEN 
    1034890                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 
     
    1053909                  ! - Maximum V-velocity difference 
    1054910                  zverr(:,:)   = 0._wp    
    1055                   DO jj = 2, jpj - 1 
    1056                      DO ji = 2, jpi - 1 
     911                  DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     912                   
    1057913                        zverr(ji,jj) = ABS ( ( v_ice(ji,jj) - zv_b(ji,jj) ) ) * vmask(ji,jj,1) 
    1058                      END DO 
    1059                   END DO 
    1060914                   
     915                  END_2D 
     916                            
    1061917                  zverr_max = MAXVAL( zverr ) 
    1062918                  CALL mpp_max( 'icedyn_rhg_evp', zverr_max )   ! max over the global domain - damned! 
     
    1083939            ! 
    1084940            !--------------------------------------------------------------------------------------- 
    1085  
    1086             IF( nn_rhg_chkcvg/=0 .AND. MOD ( i_inn - 1, nn_vp_chkcvg ) == 0 ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 
    1087                       &                         zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 
    1088  
    1089             IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' 
     941            IF( nn_rhg_chkcvg/=0 .AND. MOD ( i_inn - 1, nn_vp_chkcvg ) == 0 ) THEN 
     942 
     943               CALL rhg_cvg_vp( kt, i_out, i_inn, i_inn_tot, nn_vp_nout, nn_vp_ninn, nn_nvp,        & 
     944                      &         u_ice, v_ice, zu_b, zv_b, zu_c, zv_c,                               & 
     945                      &         zmt, za_iU, za_iV, zuerr_max, zverr_max, zglob_area,                & 
     946                      &         zrhsu, zAU, zBU, zCU, zDU, zEU, zFU,                                & 
     947                      &         zrhsv, zAV, zBV, zCV, zDV, zEV, zFV,                                & 
     948                                zvel_res, zvel_diff ) 
     949 
     950            ENDIF 
    1090951 
    1091952         END DO ! i_inn, end of inner loop 
     
    1093954      END DO ! End of outer loop (i_out) ============================================================================================= 
    1094955 
    1095       IF ( lwp ) WRITE(numout,*) ' We are out of outer loop ' 
    1096  
    1097       CALL lbc_lnk( 'icedyn_rhg_vp', zFU  , 'U',  1., zFV  , 'V',  1. ) 
    1098       CALL lbc_lnk( 'icedyn_rhg_vp', zBU_prime  , 'U',  1., zBV_prime  , 'V',  1. ) 
    1099       CALL lbc_lnk( 'icedyn_rhg_vp', zFU_prime  , 'U',  1., zFV_prime  , 'V',  1. ) 
    1100       CALL lbc_lnk( 'icedyn_rhg_vp', zCU_prime  , 'U',  1., zCV_prime  , 'V',  1. ) 
    1101  
    1102 !!$      CALL iom_put( 'zFU'           , zFU            ) ! MV DEBUG 
    1103 !!$      CALL iom_put( 'zBU_prime'     , zBU_prime      ) ! MV DEBUG 
    1104 !!$      CALL iom_put( 'zCU_prime'     , zCU_prime      ) ! MV DEBUG 
    1105 !!$      CALL iom_put( 'zFU_prime'     , zFU_prime      ) ! MV DEBUG 
    1106 !!$ 
    1107 !!$      CALL iom_put( 'zFV'           , zFV            ) ! MV DEBUG 
    1108 !!$      CALL iom_put( 'zBV_prime'     , zBV_prime      ) ! MV DEBUG 
    1109 !!$      CALL iom_put( 'zCV_prime'     , zCV_prime      ) ! MV DEBUG 
    1110 !!$      CALL iom_put( 'zFV_prime'     , zFV_prime      ) ! MV DEBUG 
    1111  
    1112       CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
    1113  
    1114 !!$      IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg ' 
    1115 !!$      IF( iom_use('uice_dbg' ) )   CALL iom_put( 'uice_dbg'   , u_ice    )                              ! ice velocity u after solver 
    1116 !!$      IF( iom_use('vice_dbg' ) )   CALL iom_put( 'vice_dbg'   , v_ice    )                              ! ice velocity v after solver 
    1117              
     956      IF( nn_rhg_chkcvg/=0  ) THEN 
     957           
     958         IF( iom_use('velo_res') )   CALL iom_put( 'velo_res', zvel_res  )   ! linear system residual  @last inner&outer iteration 
     959         IF( iom_use('velo_ero') )   CALL iom_put( 'velo_ero', zvel_diff )   ! abs velocity difference @last outer iteration 
     960         IF( iom_use('uice_eri') )   CALL iom_put( 'uice_eri', zuerr     )   ! abs velocity difference @last inner iteration 
     961         IF( iom_use('vice_eri') )   CALL iom_put( 'vice_eri', zverr     )   ! abs velocity difference @last inner iteration 
     962 
     963         DEALLOCATE( zvel_res , zvel_diff ) 
     964         
     965      ENDIF ! nn_rhg_chkcvg 
     966 
    1118967      !------------------------------------------------------------------------------! 
    1119968      ! 
    1120       ! --- Convergence diagnostics  
     969      ! --- Recompute delta, shear and div (inputs for mechanical redistribution)  
    1121970      ! 
    1122971      !------------------------------------------------------------------------------! 
    1123  
    1124       IF( nn_rhg_chkcvg /= 0 ) THEN 
    1125            
    1126          IF( iom_use('uice_cvg')  ) THEN 
    1127             CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , &                ! ice velocity difference at last iteration 
    1128                   &                        ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
    1129          ENDIF    
    1130          
    1131       ENDIF 
    1132  
    1133 !!$      ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead 
    1134 !!$      DO jj = 2, jpj - 1 
    1135 !!$         DO ji = 2, jpi - 1    
    1136 !!$ 
    1137 !!$             u_ice(ji,jj) =   zmsk00x(ji,jj)                               &  
    1138 !!$      &         * (           zmsk01x(ji,jj)   * u_oce(ji,jj) * 0.01_wp    & 
    1139 !!$                  + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) 
    1140 !!$ 
    1141 !!$             v_ice(ji,jj) =   zmsk00y(ji,jj)                               &  
    1142 !!$      &         * (           zmsk01y(ji,jj)   * v_oce(ji,jj) * 0.01_wp    & 
    1143 !!$                  + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) 
    1144 !!$ 
    1145 !!$         END DO 
    1146 !!$      END DO 
    1147 !!$ 
    1148 !!$      CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
    1149 !!$ 
    1150 !!$      IF ( lwp ) WRITE(numout,*) ' Velocity replaced ' 
    1151  
    1152       ! END DEBUG 
    1153  
    1154       !------------------------------------------------------------------------------! 
    1155       ! 
    1156       ! --- Recompute delta, shear and div (inputs for mechanical redistribution)  
    1157       ! 
    1158       !------------------------------------------------------------------------------! 
    1159972      ! 
    1160973      ! MV OPT: subroutinize ? 
    1161  
    1162       DO jj = 1, jpj - 1 
    1163          DO ji = 1, jpi - 1 
     974      DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) ! 1->jpj-1; 1->jpi-1 
    1164975 
    1165976            ! shear at F points 
     
    1168979               &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj) 
    1169980 
    1170          END DO 
    1171       END DO            
    1172        
    1173       DO jj = 2, jpj - 1 
    1174          DO ji = 2, jpi - 1 !  
     981      END_2D       
     982       
     983      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
    1175984             
    1176985            ! tension**2 at T points 
     
    1180989            zdt2 = zdt * zdt 
    1181990             
    1182             zten_i(ji,jj) = zdt 
     991            ztens(ji,jj)    = zdt 
    1183992             
    1184993            ! shear**2 at T points (doc eq. A16) 
     
    1187996               &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
    1188997             
    1189             ! shear at T points 
    1190             pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     998            ! maximum shear rate at T points (includees tension, output only) 
     999            pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) ! i think this is maximum shear rate and not actual shear --- i'm not totally sure here 
     1000 
     1001            ! shear at T-points 
     1002            zshear(ji,jj)   = SQRT( zds2 ) 
    11911003 
    11921004            ! divergence at T points 
     
    11941006               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    11951007               &             ) * r1_e1e2t(ji,jj) 
     1008 
     1009            zdiv2          =  pdivu_i(ji,jj) *  pdivu_i(ji,jj) 
    11961010             
    11971011            ! delta at T points 
    1198             zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
     1012            zdelta         = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
    11991013            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    1200              
    1201             !pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
    1202             pdelta_i(ji,jj) = zdelta + rn_creepl 
    1203  
    1204          END DO 
    1205       END DO 
    1206  
    1207       IF ( lwp ) WRITE(numout,*) ' Deformation recalculated ' 
     1014 
     1015            pdelta_i(ji,jj) = zdelta + rn_creepl ! * rswitch 
     1016 
     1017      END_2D 
    12081018       
    12091019      CALL lbc_lnk( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
     
    12221032      ! 
    12231033      ! ---- Sea ice stresses at T-points 
    1224       IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
    1225        
    1226          DO jj = 2, jpj - 1 
    1227             DO ji = 2, jpi - 1 
    1228                zp_deltastar_t(ji,jj)   =   strength(ji,jj) / pdelta_i(ji,jj)  
    1229                zfac                    =   zp_deltastar_t(ji,jj)  
     1034      IF ( iom_use('normstr')    .OR. iom_use('sheastr')    .OR. & 
     1035     &     iom_use('intstrx')    .OR. iom_use('intstry')    .OR. & 
     1036     &     iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
     1037       
     1038         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1039          
     1040               zp_delstar_t(ji,jj)     =   strength(ji,jj) / pdelta_i(ji,jj)  
     1041               zfac                    =   zp_delstar_t(ji,jj)  
    12301042               zs1(ji,jj)              =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 
    1231                zs2(ji,jj)              =   zfac * z1_ecc2 * zten_i(ji,jj) 
    1232                zs12(ji,jj)             =   zfac * z1_ecc2 * pshear_i(ji,jj) 
    1233             END DO 
    1234          END DO 
     1043               zs2(ji,jj)              =   zfac * z1_ecc2 * ztens(ji,jj) 
     1044               zs12(ji,jj)             =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bug 12 nov 
     1045 
     1046         END_2D 
    12351047 
    12361048         CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. ) 
     
    12411053      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
    12421054 
    1243          DO jj = 1, jpj - 1 
    1244             DO ji = 1, jpi - 1 
     1055         DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) ! 1->jpj-1; 1->jpi-1 
    12451056             
    12461057               ! P/delta* at F points 
    1247                zp_deltastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) ) 
     1058               zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) ) 
    12481059                
    12491060               ! s12 at F-points  
    1250                zs12f(ji,jj) = zp_deltastar_f * z1_ecc2 * zds(ji,jj) 
     1061               zs12f(ji,jj) = zp_delstar_f * z1_ecc2 * zds(ji,jj) 
    12511062                
    1252             END DO 
    1253          END DO 
     1063         END_2D 
    12541064 
    12551065         CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. ) 
    12561066          
    12571067      ENDIF 
    1258  
    1259       IF ( lwp ) WRITE(numout,*) ' zs12f recalculated ' 
    12601068 
    12611069      ! 
     
    12711079 
    12721080         !--- Recalculate oceanic stress at last inner iteration 
    1273          DO jj = 2, jpj - 1 
    1274             DO ji = 2, jpi - 1 
     1081         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
    12751082 
    12761083                !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 
     
    12881095                ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
    12891096                 
    1290             END DO 
    1291          END DO 
     1097         END_2D 
    12921098          
    12931099         ! 
    12941100         CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, & 
    1295 !            &                           ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 
     1101!            &                          ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. ) 
    12961102         ! 
    12971103         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 ) 
     
    13081114      ! --- Divergence, shear and strength --- ! 
    13091115      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
    1310       IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
     1116      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! maximum shear rate 
    13111117      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta 
    13121118      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    13131119 
    1314       IF ( lwp ) WRITE(numout,*) 'Some terms recalculated ' 
    1315  
    13161120      ! --- Stress tensor invariants (SIMIP diags) --- ! 
    13171121      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
     
    13251129         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    13261130         !          
    1327          DO jj = 2, jpj - 1 
    1328             DO ji = 2, jpi - 1 
     1131         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
    13291132               ! Stress invariants 
    1330                zsig_I(ji,jj)    =   zs1(ji,jj) * 0.5_wp                                        ! 1st invariant, aka average normal stress aka negative pressure 
    1331                zsig_II(ji,jj)   =   SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj) )  ! 2nd invariant, aka maximum shear stress                
    1332             END DO 
    1333          END DO 
     1133               zsig_I(ji,jj)    =   zs1(ji,jj) * 0.5_wp   ! 1st invariant, aka average normal stress aka negative pressure 
     1134               zsig_II(ji,jj)   =   0.5_wp * SQRT ( zs2(ji,jj) * zs2(ji,jj) + 4. * zs12(ji,jj) * zs12(ji,jj) )   ! 2nd invariant, aka maximum shear stress                
     1135         END_2D 
    13341136 
    13351137         CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.) 
     
    13411143          
    13421144      ENDIF 
    1343  
    1344       IF ( lwp ) WRITE(numout,*) 'SIMIP work done' 
    13451145 
    13461146      ! --- Normalized stress tensor principal components --- ! 
     
    13551155         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    13561156         !          
    1357          DO jj = 2, jpj - 1 
    1358             DO ji = 2, jpi - 1 
     1157         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1158          
    13591159               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates  
    13601160               !                        and **deformations** at current iterates 
    13611161               !                        following Lemieux & Dupont (2020) 
    1362                zfac             =   zp_deltastar_t(ji,jj) 
    1363                zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) ) 
    1364 !!$               zsig1            = 0._wp !!! FUCKING DEBUG TEST !!! 
    1365                zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj) 
    1366                zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj) 
     1162               zfac             =   zp_delstar_t(ji,jj) 
     1163               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltat(ji,jj) ) 
     1164               zsig2            =   zfac * z1_ecc2 * ztens(ji,jj) 
     1165               zsig12           =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bugfix 12 Nov 
    13671166                
    13681167               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
    1369                zsig_I(ji,jj)    =   zsig1 * 0.5_wp                              ! 1st invariant 
    1370                zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd invariant 
     1168               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                       ! 1st invariant 
     1169               zsig_II(ji,jj)   =   0.5_wp * SQRT ( zsig2 * zsig2 + 4. *zsig12 * zsig12 )   ! 2nd invariant 
    13711170 
    13721171               ! Normalized  principal stresses (used to display the ellipse) 
     
    13741173               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
    13751174               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
    1376             END DO 
    1377          END DO 
    1378          IF ( lwp ) WRITE(numout,*) 'Some shitty stress work done' 
     1175                
     1176         END_2D 
    13791177         ! 
    13801178         CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.) 
    1381          !       
    1382          IF ( lwp ) WRITE(numout,*) ' Beauaaaarflblbllll ' 
    13831179         ! 
    13841180         CALL iom_put( 'sig1_pnorm' , zsig1_p )  
     
    13861182 
    13871183         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II ) 
    1388  
    1389          IF ( lwp ) WRITE(numout,*) ' So what ??? ' 
    13901184          
    13911185      ENDIF 
     
    13961190 
    13971191         ! --- Recalculate Coriolis stress at last inner iteration 
    1398          DO jj = 2, jpj - 1 
    1399             DO ji = 2, jpi - 1 
     1192         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
    14001193                ! --- U-component  
    14011194                zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  & 
     
    14051198                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    14061199                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    1407             END DO 
    1408          END DO 
     1200         END_2D 
    14091201         ! 
    14101202         CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., & 
     
    14211213 
    14221214         ! Recalculate internal forces (divergence of stress tensor) at last inner iteration 
    1423          DO jj = 2, jpj - 1 
    1424             DO ji = 2, jpi - 1 
     1215         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1216 
    14251217               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
    14261218                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     
    14291221                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
    14301222                  &                  ) * r1_e1e2u(ji,jj) 
     1223 
    14311224               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
    14321225                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     
    14351228                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
    14361229                  &                  ) * r1_e1e2v(ji,jj) 
    1437             END DO 
    1438          END DO 
     1230 
     1231         END_2D 
    14391232             
    14401233         CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. ) 
     
    14521245            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 
    14531246         ! 
    1454          DO jj = 2, jpj - 1 
    1455             DO ji = 2, jpi - 1 
    1456                ! 2D ice mass, snow mass, area transport arrays (X, Y) 
     1247         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1               ! 2D ice mass, snow mass, area transport arrays (X, Y) 
     1248 
    14571249               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 
    14581250               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj) 
     
    14671259               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   '' 
    14681260 
    1469             END DO 
    1470          END DO 
    1471  
     1261         END_2D 
     1262          
    14721263         CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., & 
    14731264            &                           zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., & 
     
    14891280    
    14901281    
    1491     
    1492    SUBROUTINE rhg_cvg_vp( kt, kiter, kitermax, pu, pv, pmt, puerr_max, pverr_max, pglob_area, & 
    1493                   &       prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV ) 
    1494     
     1282   SUBROUTINE rhg_cvg_vp( kt, kitout, kitinn, kitinntot, kitoutmax, kitinnmax, kitinntotmax , & 
     1283                  &       pu, pv, pub, pvb, pub_outer, pvb_outer                     , & 
     1284                  &       pmt, pat_iu, pat_iv, puerr_max, pverr_max, pglob_area      , & 
     1285                  &       prhsu, pAU, pBU, pCU, pDU, pEU, pFU                        , & 
     1286                  &       prhsv, pAV, pBV, pCV, pDV, pEV, pFV                        , &    
     1287                  &       pvel_res, pvel_diff                                            ) 
     1288      !! 
    14951289      !!---------------------------------------------------------------------- 
    14961290      !!                    ***  ROUTINE rhg_cvg_vp  *** 
     
    15071301      !! 
    15081302      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)    
     1303      !! 
    15091304      !!---------------------------------------------------------------------- 
    1510       INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax      ! ocean time-step index 
    1511       REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pmt              ! now velocity and mass per unit area  
    1512       REAL(wp),                 INTENT(in) ::   puerr_max, pverr_max     ! absolute mean velocity difference 
    1513       REAL(wp),                 INTENT(in) ::   pglob_area               ! global ice area 
    1514       REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsu, pAU, pBU, pCU, pDU, pEU ! linear system coefficients  
    1515       REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsv, pAV, pBV, pCV, pDV, pEV 
    1516       !! 
    1517       INTEGER           ::   it, idtime, istatus, ix_dim, iy_dim 
     1305      !! 
     1306      INTEGER ,                 INTENT(in) ::   kt, kitout, kitinn, kitinntot    ! ocean model iterate, outer, inner and total n-iterations 
     1307      INTEGER ,                 INTENT(in) ::   kitoutmax, kitinnmax             ! max number of outer & inner iterations 
     1308      INTEGER ,                 INTENT(in) ::   kitinntotmax                     ! max number of total sub-iterations 
     1309      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb                 ! now & sub-iter-before velocities 
     1310      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pub_outer, pvb_outer             ! velocities @before outer iterations 
     1311      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pmt, pat_iu, pat_iv              ! mass at T-point, ice concentration at U&V 
     1312      REAL(wp),                 INTENT(in) ::   puerr_max, pverr_max             ! absolute mean velocity difference 
     1313      REAL(wp),                 INTENT(in) ::   pglob_area                       ! global ice area 
     1314      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsu, pAU, pBU, pCU, pDU, pEU, pFU ! linear system coefficients  
     1315      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsv, pAV, pBV, pCV, pDV, pEV, pFV 
     1316      REAL(wp), DIMENSION(:,:), INTENT(inout) ::  pvel_res                       ! velocity residual @last inner iteration 
     1317      REAL(wp), DIMENSION(:,:), INTENT(inout) ::  pvel_diff                      ! velocity difference @last outer iteration 
     1318      !! 
     1319 
     1320      INTEGER           ::   idtime, istatus, ix_dim, iy_dim 
    15181321      INTEGER           ::   ji, jj          ! dummy loop indices 
    1519       REAL(wp)          ::   zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars 
    1520       REAL(wp)          ::   z1_pglob_area 
     1322      INTEGER           ::   it_inn_file, it_out_file 
     1323      REAL(wp)          ::   zu_res_mean, zv_res_mean, zvel_res_mean                  ! mean residuals of the linear system 
     1324      REAL(wp)          ::   zu_mad, zv_mad, zvel_mad                                 ! mean absolute deviation, sub-iterates 
     1325      REAL(wp)          ::   zu_mad_outer, zv_mad_outer, zvel_mad_outer               ! mean absolute deviation, outer-iterates 
     1326      REAL(wp)          ::   zvel_err_max, zmke, zu, zv                               ! local scalars 
     1327      REAL(wp)          ::   z1_pglob_area                                            ! inverse global ice area 
     1328 
    15211329      REAL(wp), DIMENSION(jpi,jpj) ::   zu_res, zv_res, zvel2                         ! local arrays 
     1330      REAL(wp), DIMENSION(jpi,jpj) ::   zu_diff, zv_diff                              ! local arrays 
    15221331                                                                              
    15231332      CHARACTER(len=20) ::   clname 
    15241333      !!---------------------------------------------------------------------- 
    15251334 
     1335 
    15261336      IF( lwp ) THEN 
     1337 
    15271338         WRITE(numout,*) 
    15281339         WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 
    15291340         WRITE(numout,*) '~~~~~~~~~~~' 
    1530          WRITE(numout,*) ' kiter    =  : ', kiter 
    1531          WRITE(numout,*) ' kitermax =  : ', kitermax 
     1341         WRITE(numout,*) ' kt          =  : ', kt 
     1342         WRITE(numout,*) ' kitout      =  : ', kitout 
     1343         WRITE(numout,*) ' kitinn      =  : ', kitinn 
     1344         WRITE(numout,*) ' kitinntot   =  : ', kitinntot 
     1345         WRITE(numout,*) ' kitoutmax (nn_vp_nout) =  ', kitoutmax 
     1346         WRITE(numout,*) ' kitinnmax (nn_vp_ninn) =  ', kitinnmax 
     1347         WRITE(numout,*) ' kitinntotmax (nn_nvp)  =  ', kitinntotmax 
     1348         WRITE(numout,*) 
     1349 
    15321350      ENDIF 
    15331351 
     1352      z1_pglob_area = 1._wp / pglob_area      ! inverse global ice area 
     1353 
    15341354      ! create file 
    1535       IF( kt == nit000 .AND. kiter == 1 ) THEN 
     1355      IF( kt == nit000 .AND. kitinntot == 1 ) THEN 
    15361356         ! 
    15371357         IF( lwm ) THEN 
     
    15451365            istatus = NF90_DEF_DIM( ncvgid, 'y'     , jpj, iy_dim ) 
    15461366 
    1547             ! i suggest vel_dif instead 
    1548             istatus = NF90_DEF_VAR( ncvgid, 'u_res'  , NF90_DOUBLE  , (/ idtime /), nvarid_ures ) 
    1549             istatus = NF90_DEF_VAR( ncvgid, 'v_res'  , NF90_DOUBLE  , (/ idtime /), nvarid_vres ) 
    1550             istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE  , (/ idtime /), nvarid_velres ) 
    1551             istatus = NF90_DEF_VAR( ncvgid, 'u_dif'  , NF90_DOUBLE  , (/ idtime /), nvarid_udif ) 
    1552             istatus = NF90_DEF_VAR( ncvgid, 'v_dif'  , NF90_DOUBLE  , (/ idtime /), nvarid_vdif ) 
    1553             istatus = NF90_DEF_VAR( ncvgid, 'vel_dif', NF90_DOUBLE  , (/ idtime /), nvarid_veldif ) 
     1367            istatus = NF90_DEF_VAR( ncvgid, 'u_res'         , NF90_DOUBLE  , (/ idtime /), nvarid_ures ) 
     1368            istatus = NF90_DEF_VAR( ncvgid, 'v_res'         , NF90_DOUBLE  , (/ idtime /), nvarid_vres ) 
     1369            istatus = NF90_DEF_VAR( ncvgid, 'vel_res'       , NF90_DOUBLE  , (/ idtime /), nvarid_velres ) 
     1370 
     1371            istatus = NF90_DEF_VAR( ncvgid, 'uerr_max_sub'  , NF90_DOUBLE  , (/ idtime /), nvarid_uerr_max ) 
     1372            istatus = NF90_DEF_VAR( ncvgid, 'verr_max_sub'  , NF90_DOUBLE  , (/ idtime /), nvarid_verr_max ) 
     1373            istatus = NF90_DEF_VAR( ncvgid, 'velerr_max_sub', NF90_DOUBLE  , (/ idtime /), nvarid_velerr_max ) 
     1374 
     1375            istatus = NF90_DEF_VAR( ncvgid, 'umad_sub'      , NF90_DOUBLE  , (/ idtime /), nvarid_umad ) 
     1376            istatus = NF90_DEF_VAR( ncvgid, 'vmad_sub'      , NF90_DOUBLE  , (/ idtime /), nvarid_vmad ) 
     1377            istatus = NF90_DEF_VAR( ncvgid, 'velmad_sub'    , NF90_DOUBLE  , (/ idtime /), nvarid_velmad ) 
     1378             
     1379            istatus = NF90_DEF_VAR( ncvgid, 'umad_outer'    , NF90_DOUBLE  , (/ idtime /), nvarid_umad_outer   ) 
     1380            istatus = NF90_DEF_VAR( ncvgid, 'vmad_outer'    , NF90_DOUBLE  , (/ idtime /), nvarid_vmad_outer   ) 
     1381            istatus = NF90_DEF_VAR( ncvgid, 'velmad_outer'  , NF90_DOUBLE  , (/ idtime /), nvarid_velmad_outer ) 
     1382 
    15541383            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE  , (/ idtime /), nvarid_mke ) 
    15551384 
    1556             istatus = NF90_DEF_VAR( ncvgid, 'u_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_ures_xy) 
    1557             istatus = NF90_DEF_VAR( ncvgid, 'v_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_vres_xy) 
    1558  
    15591385            istatus = NF90_ENDDEF(ncvgid) 
    15601386 
     
    15631389      ENDIF 
    15641390 
    1565       IF ( lwp ) WRITE(numout,*) ' File created ' 
    1566  
    1567       ! --- Max absolute velocity difference with previous iterate (zveldif) 
    1568       zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code  
    1569                                             ! if puerrmask and pverrmax are masked at 15% (TEST) 
    1570  
    1571       ! ---  Mean residual and kinetic energy 
    1572       IF ( kiter == 1 ) THEN 
    1573  
    1574          zu_res_mean = 0._wp 
    1575          zv_res_mean = 0._wp 
    1576          zvelres     = 0._wp 
    1577          zmke        = 0._wp 
     1391      !------------------------------------------------------------ 
     1392      ! 
     1393      ! Max absolute velocity difference with previous sub-iterate 
     1394      ! ( zvel_err_max ) 
     1395      ! 
     1396      !------------------------------------------------------------ 
     1397      ! 
     1398      ! This comes from the criterion used to stop the iterative procedure 
     1399      zvel_err_max   = 0.5_wp * ( puerr_max + pverr_max ) ! average of U- and V- maximum error over the whole domain 
     1400 
     1401      !---------------------------------------------- 
     1402      ! 
     1403      ! Mean-absolute-deviation (sub-iterates) 
     1404      ! ( zu_mad, zv_mad, zvel_mad) 
     1405      ! 
     1406      !---------------------------------------------- 
     1407      ! 
     1408      ! U 
     1409      zu_diff(:,:) = 0._wp 
     1410       
     1411      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1412       
     1413         zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * z1_pglob_area 
     1414       
     1415      END_2D 
     1416       
     1417      ! V 
     1418      zv_diff(:,:)   = 0._wp 
     1419       
     1420      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1421       
     1422         zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * z1_pglob_area 
     1423       
     1424      END_2D 
     1425 
     1426      ! global sum & U-V average 
     1427      CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_diff,  'U',  1., zv_diff , 'V',  1. ) 
     1428      zu_mad   = glob_sum( 'icedyn_rhg_vp : ', zu_diff ) 
     1429      zv_mad   = glob_sum( 'icedyn_rhg_vp : ', zv_diff ) 
     1430 
     1431      zvel_mad = 0.5_wp * ( zu_mad + zv_mad ) 
     1432 
     1433      !----------------------------------------------- 
     1434      ! 
     1435      ! Mean-absolute-deviation (outer-iterates) 
     1436      ! ( zu_mad_outer, zv_mad_outer, zvel_mad_outer) 
     1437      ! 
     1438      !----------------------------------------------- 
     1439      ! 
     1440      IF ( kitinn == kitinnmax ) THEN ! only work at the end of outer iterates  
     1441 
     1442         ! * U 
     1443         zu_diff(:,:) = 0._wp 
     1444          
     1445         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1446          
     1447            zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub_outer(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * & 
     1448                              &    z1_pglob_area 
     1449                               
     1450         END_2D 
     1451          
     1452         ! * V 
     1453         zv_diff(:,:)   = 0._wp 
     1454          
     1455         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1456          
     1457            zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb_outer(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * & 
     1458                              &    z1_pglob_area 
     1459          
     1460         END_2D 
     1461          
     1462         ! Global ice-concentration, grid-cell-area weighted mean 
     1463         CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_diff,  'U',  1., zv_diff , 'V',  1. ) ! abs behaves as a scalar no ? 
     1464 
     1465         zu_mad_outer   = glob_sum( 'icedyn_rhg_vp : ', zu_diff ) 
     1466         zv_mad_outer   = glob_sum( 'icedyn_rhg_vp : ', zv_diff ) 
     1467    
     1468         ! Average of both U & V 
     1469         zvel_mad_outer = 0.5_wp * ( zu_mad_outer + zv_mad_outer ) 
     1470                   
     1471      ENDIF 
     1472 
     1473      ! --- Spatially-resolved absolute difference to send back to main routine  
     1474      ! (last iteration only, T-point) 
     1475 
     1476      IF ( kitinntot == kitinntotmax) THEN 
     1477 
     1478         zu_diff(:,:) = 0._wp 
     1479         zv_diff(:,:) = 0._wp 
     1480 
     1481         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1482 
     1483               zu_diff(ji,jj) = ( ABS ( ( pu(ji-1,jj) - pub_outer(ji-1,jj) ) ) * umask(ji-1,jj,1) & 
     1484    &                           + ABS ( ( pu(ji,jj  ) - pub_outer(ji,jj)   ) ) * umask(ji,jj,1) ) & 
     1485    &                           / ( umask(ji-1,jj,1) + umask(ji,jj,1) ) 
     1486 
     1487               zv_diff(ji,jj) = ( ABS ( ( pv(ji,jj-1)   - pvb_outer(ji,jj-1) ) ) * vmask(ji,jj-1,1) & 
     1488    &                           + ABS ( ( pv(ji,jj  ) - pvb_outer(ji,jj)     ) ) * vmask(ji,jj,1)   & 
     1489    &                           / ( vmask(ji,jj-1,1) + vmask(ji,jj,1) ) ) 
     1490      
     1491    
     1492         END_2D 
     1493          
     1494         CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_diff,  'T',  1., zv_diff , 'T',  1. ) 
     1495         pvel_diff(:,:) = 0.5_wp * ( zu_diff(:,:) + zv_diff(:,:) ) 
    15781496 
    15791497      ELSE 
    15801498 
    1581       ! -- Mean residual (N/m^2), zu_res_mean 
    1582       ! Here we take the residual of the linear system (N/m^2),  
    1583       ! We define it as in mitgcm: square-root of area-weighted mean square residual 
    1584       ! Local residual r = Ax - B expresses to which extent the momentum balance is verified  
    1585       ! i.e., how close we are to a solution 
    1586  
    1587       IF ( lwp ) WRITE(numout,*) ' TEST 1 '  
    1588  
    1589       z1_pglob_area = 1._wp / pglob_area 
    1590  
    1591       zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 
    1592  
    1593       DO jj = 2, jpj - 1 
    1594          DO ji = 2, jpi - 1                                       
    1595             zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               & 
    1596                &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 
    1597                             
    1598             zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               & 
    1599                &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 
    1600  
    1601             zu_res(ji,jj)  = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * e1e2u(ji,jj) * z1_pglob_area 
    1602             zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area 
    1603  
    1604          END DO 
    1605       END DO                   
    1606  
    1607       IF ( lwp ) WRITE(numout,*) ' TEST 2 '  
    1608       zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 
    1609       zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 
    1610       IF ( lwp ) WRITE(numout,*) ' TEST 3 '  
    1611       zvelres     = 0.5_wp * ( zu_res_mean + zv_res_mean ) 
    1612  
    1613       IF ( lwp ) WRITE(numout,*) ' TEST 4 '  
    1614                                            
    1615       ! -- Global mean kinetic energy per unit area (J/m2) 
    1616       zvel2(:,:) = 0._wp 
    1617       DO jj = 2, jpj - 1 
    1618          DO ji = 2, jpi - 1                    
    1619             zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 
    1620             zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 
    1621             zvel2(ji,jj)  = zu*zu + zv*zv              ! square of ice velocity at T-point   
    1622          END DO 
    1623       END DO 
    1624         
    1625       IF ( lwp ) WRITE(numout,*) ' TEST 5 '  
    1626  
    1627       zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 
    1628  
    1629       IF ( lwp ) WRITE(numout,*) ' TEST 6 ' 
    1630  
    1631       ENDIF ! kiter 
     1499         pvel_diff(:,:) = 0._wp 
     1500 
     1501      ENDIF 
     1502 
     1503      !--------------------------------------- 
     1504      ! 
     1505      ! ---  Mean residual & kinetic energy 
     1506      ! 
     1507      !--------------------------------------- 
     1508 
     1509      IF ( kitinntot == 1 ) THEN 
     1510 
     1511         zu_res_mean   = 0._wp 
     1512         zv_res_mean   = 0._wp 
     1513         zvel_res_mean = 0._wp 
     1514         zmke          = 0._wp 
     1515 
     1516      ELSE 
     1517 
     1518         ! * Mean residual (N/m2) 
     1519         ! Here we take the residual of the linear system (N/m2),  
     1520         ! We define it as in mitgcm: global area-weighted mean of square-root residual 
     1521         ! Local residual r = Ax - B expresses to which extent the momentum balance is verified  
     1522         ! i.e., how close we are to a solution 
     1523         zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 
     1524          
     1525         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1526 
     1527               zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               & 
     1528                  &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 
     1529               zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               & 
     1530                  &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 
     1531 
     1532!              zu_res(ji,jj)  = pFU(ji,jj) - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) 
     1533!              zv_res(ji,jj)  = pFV(ji,jj) - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) 
     1534    
     1535               zu_res(ji,jj)  = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * pat_iu(ji,jj) * e1e2u(ji,jj) * z1_pglob_area 
     1536               zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * pat_iv(ji,jj) * e1e2v(ji,jj) * z1_pglob_area 
     1537    
     1538         END_2D 
     1539          
     1540         ! Global ice-concentration, grid-cell-area weighted mean 
     1541         CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_res,  'U', 1., zv_res , 'V', 1. ) 
     1542    
     1543         zu_res_mean   = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 
     1544         zv_res_mean   = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 
     1545         zvel_res_mean = 0.5_wp * ( zu_res_mean + zv_res_mean ) 
     1546    
     1547         ! --- Global mean kinetic energy per unit area (J/m2) 
     1548         zvel2(:,:) = 0._wp 
     1549 
     1550         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
    16321551                   
    1633          !                                                ! ==================== ! 
    1634  
    1635       ! time 
    1636       it = ( kt - nit000 ) * kitermax + kiter 
    1637  
    1638  
     1552               zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 
     1553               zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 
     1554               zvel2(ji,jj)  = zu*zu + zv*zv              ! square of ice velocity at T-point   
     1555 
     1556         END_2D 
     1557                    
     1558         zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 
     1559    
     1560      ENDIF ! kitinntot 
     1561 
     1562      !--- Spatially-resolved residual at last iteration to send back to main routine (last iteration only) 
     1563      !--- Calculation @T-point 
     1564 
     1565      IF ( kitinntot == kitinntotmax) THEN 
     1566 
     1567         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1568 
     1569               zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               & 
     1570                  &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 
     1571               zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               & 
     1572                  &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 
     1573 
     1574               zu_res(ji,jj)  = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1)  
     1575               zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1)  
     1576 
     1577         END_2D 
     1578          
     1579         CALL lbc_lnk( 'icedyn_rhg_cvg_vp', zu_res,  'U',  1., zv_res , 'V',  1. ) 
     1580 
     1581         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! 2->jpj-1; 2->jpi-1 
     1582          
     1583               pvel_res(ji,jj) = 0.25_wp * ( zu_res(ji-1,jj) + zu_res(ji,jj) + zv_res(ji,jj-1) + zv_res(ji,jj) ) 
     1584          
     1585         END_2D 
     1586         CALL lbc_lnk( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. ) 
     1587 
     1588      ELSE 
     1589 
     1590         pvel_res(:,:) = 0._wp 
     1591 
     1592      ENDIF 
     1593                   
     1594      !                                                ! ==================== ! 
     1595 
     1596      it_inn_file =  ( kt - nit000 ) * kitinntotmax + kitinntot ! time step in the file 
     1597      it_out_file =  ( kt - nit000 ) * kitoutmax    + kitout 
     1598 
     1599      ! write variables 
    16391600      IF( lwm ) THEN 
    1640          ! write variables 
    1641          istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 
    1642          istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 
    1643          istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) )   ! average of u- and v- residuals 
    1644          istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) )   ! max U velocity difference, inner iterations 
    1645          istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) )   ! max V velocity difference, inner iterations 
    1646          istatus = NF90_PUT_VAR( ncvgid, nvarid_veldif, (/zveldif/), (/it/), (/1/) )   ! max U or V velocity diff between subiterations 
    1647          istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) )         ! mean kinetic energy 
    1648  
    1649          ! 
    1650          IF ( kiter == kitermax ) THEN 
    1651             WRITE(numout,*) ' Should plot the spatially dependent residual ' 
    1652             istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) )          ! U-residual, spatially dependent 
    1653             istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) )          ! V-residual, spatially dependent 
     1601 
     1602         istatus = NF90_PUT_VAR( ncvgid, nvarid_ures  , (/zu_res_mean/), (/it_inn_file/), (/1/) )        ! Residuals of the linear system, area weighted mean 
     1603         istatus = NF90_PUT_VAR( ncvgid, nvarid_vres  , (/zv_res_mean/), (/it_inn_file/), (/1/) )        ! 
     1604         istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvel_res_mean/), (/it_inn_file/), (/1/) )      ! 
     1605 
     1606         istatus = NF90_PUT_VAR( ncvgid, nvarid_uerr_max  , (/puerr_max/), (/it_inn_file/), (/1/) )      ! Max velocit_inn_filey error, sub-it_inn_fileerates 
     1607         istatus = NF90_PUT_VAR( ncvgid, nvarid_verr_max  , (/pverr_max/), (/it_inn_file/), (/1/) )      !  
     1608         istatus = NF90_PUT_VAR( ncvgid, nvarid_velerr_max, (/zvel_err_max/), (/it_inn_file/), (/1/) )   !  
     1609 
     1610         istatus = NF90_PUT_VAR( ncvgid, nvarid_umad    , (/zu_mad/)  , (/it_inn_file/), (/1/) )         ! velocit_inn_filey MAD, area/sic-weighted, sub-it_inn_fileerates 
     1611         istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad    , (/zv_mad/)  , (/it_inn_file/), (/1/) )         !  
     1612         istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad  , (/zvel_mad/), (/it_inn_file/), (/1/) )         !  
     1613 
     1614         istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/kitinntot/), (/1/) )                    ! mean kinetic energy 
     1615 
     1616         IF ( kitinn == kitinnmax ) THEN ! only print outer mad at the end of inner loop 
     1617 
     1618            istatus = NF90_PUT_VAR( ncvgid, nvarid_umad_outer    , (/zu_mad_outer/)  , (/it_out_file/), (/1/) )   ! velocity MAD, area/sic-weighted, outer-iterates 
     1619            istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad_outer    , (/zv_mad_outer/)  , (/it_out_file/), (/1/) )   ! 
     1620            istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad_outer  , (/zvel_mad_outer/), (/it_out_file/), (/1/) )   ! 
     1621 
    16541622         ENDIF 
    16551623 
    1656          ! close file 
    1657          IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax )   istatus = NF90_CLOSE(ncvgid) 
     1624         IF( kt == nitend - nn_fsbc + 1 .AND. kitinntot == kitinntotmax )    istatus = NF90_CLOSE( ncvgid ) 
    16581625      ENDIF 
    16591626       
Note: See TracChangeset for help on using the changeset viewer.