Changeset 13522


Ignore:
Timestamp:
2020-09-25T11:01:23+02:00 (4 months ago)
Author:
vancop
Message:

Close to compiling VP

Location:
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/ice.F90

    r12919 r13522  
    150150   ! 
    151151   !                                     !!** ice-rheology namelist (namdyn_rhg) ** 
     152   ! -- evp 
    152153   LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F)  
    153    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9 
    154    REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve 
     154   REAL(wp), PUBLIC ::   rn_creepl        !: creep limit !!! common with vp 
     155   REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve !!! common with vp 
    155156   INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling 
    156157   REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp)  
    157158   LOGICAL , PUBLIC ::   ln_rhg_chkcvg    !: check ice rheology convergence  
     159   ! -- vp 
     160   INTEGER , PUBLIC ::   nn_nout_vp       !: Number of outer iterations 
     161   INTEGER , PUBLIC ::   nn_ninn_vp       !: Number of inner iterations (linear system solver) 
     162   LOGICAL , PUBLIC ::   ln_smooth_vp     !: zeta viscosity smoothing (if yes, tanh function of Lemieux and Tremblay JGR 2009) 
     163   LOGICAL , PUBLIC ::   ln_zebra_vp      !: activate zebra (solve the linear system problem every odd j-band, then one every even one) 
     164   REAL(wp), PUBLIC ::   rn_relaxu_vp     !: U-relaxation factor (MV: can probably be merged with V-factor once ok) 
     165   REAL(wp), PUBLIC ::   rn_relaxv_vp     !: V-relaxation factor  
     166   REAL(wp), PUBLIC ::   rn_uerr_max_vp   !: maximum velocity error, above which a forcing error is considered and solver is stopped 
     167   REAL(wp), PUBLIC ::   rn_uerr_min_vp   !: minimum velocity error, beyond which convergence is assumed 
     168   INTEGER , PUBLIC ::   nn_cvgchk_vp     !: Number of iterations every each convergence is checked 
    158169   ! 
    159170   !                                     !!** ice-advection namelist (namdyn_adv) ** 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg.F90

    r13493 r13522  
    3737 
    3838   ! ** namelist (namrhg) ** 
    39    LOGICAL ::   ln_rhg_EVP       ! EVP rheology 
     39   LOGICAL ::   ln_rhg_EVP       ! EVP rheology ! MV: why is it declared here while all other parameters are declared in ice.F90 
    4040   LOGICAL ::   ln_rhg_VP        ! EVP rheology 
    4141   ! 
     
    120120      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    121121      !! 
    122       NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg, & 
    123     &                       ln_rhg_VP  
     122      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_chkcvg, &           !- evp 
     123    &                       ln_rhg_VP, nn_nout_vp, nn_ninn_vp, ln_smooth_vp, ln_zebra_vp, rn_relaxu_vp, rn_relaxv_vp, rn_uerr_max_vp, rn_uerr_min_vp, nn_cvgchk_vp  !-vp 
    124124      !!------------------------------------------------------------------- 
     125       
    125126      ! 
    126127      REWIND( numnam_ice_ref )         ! Namelist namdyn_rhg in reference namelist : Ice dynamics 
     
    139140         WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP    = ', ln_rhg_EVP 
    140141         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP       = ', ln_aEVP 
    141          WRITE(numout,*) '         creep limit                                       rn_creepl     = ', rn_creepl 
    142          WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc        = ', rn_ecc 
     142         WRITE(numout,*) '         creep limit                                       rn_creepl     = ', rn_creepl     ! common with VP 
     143         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc        = ', rn_ecc        ! common with VP 
    143144         WRITE(numout,*) '         number of iterations for subcycling               nn_nevp       = ', nn_nevp 
    144145         WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast     = ', rn_relast 
    145          WRITE(numout,*) '      check convergence of rheology                        ln_rhg_chkcvg = ', ln_rhg_chkcvg 
     146         WRITE(numout,*) '      check convergence of rheology                        ln_rhg_chkcvg = ', ln_rhg_chkcvg ! common with VP 
    146147         WRITE(numout,*) ''     
    147          WRITE(numout,*) '      rheology VP     (icedyn_rhg_lsr)                     ln_rhg_VP     = ', ln_rhg_VP 
     148         WRITE(numout,*) '      rheology VP     (icedyn_rhg_VP)                      ln_rhg_VP     = ', ln_rhg_VP 
     149         WRITE(numout,*) '         number of outer iterations                        nn_nout_vp    = ', nn_nout_vp 
     150         WRITE(numout,*) '         number of inner iterations                        nn_ninn_vp    = ', nn_ninn_vp 
     151         WRITE(numout,*) '         tanh viscosity smooothing                         ln_smooth_vp  = ', ln_smooth_vp 
     152         WRITE(numout,*) '         activate zebra solver                             ln_zebra_vp   = ', ln_zebra_vp 
     153         WRITE(numout,*) '         relaxation factor for u                           rn_relaxu_vp  = ', rn_relaxu_vp 
     154         WRITE(numout,*) '         relaxation factor for v                           rn_relaxv_vp  = ', rn_relaxv_vp 
     155         WRITE(numout,*) '         maximum error on velocity                         rn_uerr_max_vp = ', rn_uerr_max_vp 
     156         WRITE(numout,*) '         velocity to decide convergence                    rn_uerr_min_vp = ', rn_uerr_min_vp 
     157         WRITE(numout,*) '         iteration step for convergence check              nn_cvgchk_vp   = ', nn_cvgchk_vp 
     158                            
    148159      ENDIF 
    149160      ! 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90

    r13493 r13522  
    11! TODOLIST 
    2  
    3 ! 
    4 ! Code residual norm of the linear system 
    5 ! Verify convergence check 
    62! 
    73! Define all symbols 
    8 ! Check uses 
     4! - Do viscosity smoothing with sum (differentiable rheology) 
     5! - Re-calculate internal force diagnostic (end of the routine) 
     6! - Do proper masking of output fileds with ice mass (zmsk00 see EVP routine) 
     7! 
     8! quality control - compare code to mitGCM 
     9! quality control - comparing EVP and VP codes 
     10!  
     11! 
     12! Commit 
    913! Compile 
    1014! 
    11 ! Check 
    12 ! --- lbclnks --> what is the policy ? 
    13 ! --- boundary conditions --> how to enforce them - is the fmask strategy sufficient ? 
    14 ! --- change of independent term in the U-V solvers ? 
    15 ! --- loop ji then jj in the V-solver, is it well-written 
     15! Clarify 
     16! --- Boundary conditions --> how to enforce them - is the fmask strategy sufficient ? 
     17! --- Swap of independent term in the U-V solvers ? 
     18! --- Is stress tensor for restart needed ? 
     19! --- Is stress tensor calculated at the end of the time step 
    1620! 
    17 ! Subroutinize tridiagonal systems ? 
     21! Test 
     22! --- Can we add the 15% mask in the convergence criteria ? 
     23! --- Try ADI for u-v solver 
     24! 
     25! Add 
     26! - Charge ellipse as an output (good one from Lemieux and Dupont) 
     27! - Bottom drag 
     28! - Tensile strength 
     29! - agrif 
     30! - bdy 
     31! 
    1832! Write documentation 
    1933! 
    20  
    21 ! Convergence quality criteria 
    22 ! - average KE (for outer loop iterations) (see Lemieux et al, 2009) 
    23 ! - residual of the linear system (for inner loop iterations) 
    24 ! - velocity difference ? 
    25  
    26  
    27 ! add bottom drag 
    28 ! add tensile strength 
    29  
    30 ! quality control - compare to mitGCM 
    31 ! quality control - comparing EVP and VP 
     34! Optimize 
     35! - subroutinize common parts (diagnostics) 
     36! - namelist: rationalize common parameters EVP/VP 
     37 
    3238 
    3339MODULE icedyn_rhg_vp 
     
    7581   PUBLIC   rhg_vp_rst       ! called by icedyn_rhg.F90 
    7682 
    77    !! * Substitutions 
    78 #  include "vectopt_loop_substitute.h90" 
    79  
    8083   !! for convergence tests 
    8184   INTEGER ::   ncvgid   ! netcdf file id 
    82    INTEGER ::   nvarid   ! netcdf variable id 
     85   INTEGER ::   nvarid_ucvg   ! netcdf variable id 
     86   INTEGER ::   nvarid_ures 
     87   INTEGER ::   nvarid_vres 
     88   INTEGER ::   nvarid_velres 
     89   INTEGER ::   nvarid_udif 
     90   INTEGER ::   nvarid_vdif 
     91   INTEGER ::   nvarid_mke 
    8392   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
    8493   !!---------------------------------------------------------------------- 
     
    8796   !! Software governed by the CeCILL license (see ./LICENSE) 
    8897   !!---------------------------------------------------------------------- 
     98    
    8999CONTAINS 
    90100 
     
    95105      !!                             VP-LSR-C-grid 
    96106      !! 
    97       !! ** purpose : determines sea ice drift from wind stress, ice-ocean 
     107      !! ** Purpose : determines sea ice drift from wind stress, ice-ocean 
    98108      !!  stress and sea-surface slope. Internal forces assume viscous-plastic rheology (Hibler, 1979) 
     109      !!  
     110      !! ** Method 
    99111      !!   
    100112      !!  The resolution algorithm  follows from Zhang and Hibler (1997) and Losch (2010) 
     113      !!  with elements from Lemieux and Tremblay (2008) and Lemieux and Tremblay (2009) 
    101114      !!   
    102       !!  It is based on a linear slope relaxation algorithm (LSR) 
    103       !! 
    104       !!  * Time stepping 
    105       !! 
    106       !!  The idea is to arrange the momentum equations as follows: 
     115      !!  The components of the momentum equations are arranged following the ideas of Zhang and Hibler (1997) 
     116      !! 
    107117      !!  f1(u) = g1(v) 
    108118      !!  f2(v) = g2(v) 
    109119      !! 
    110       !!  and use two sub-time-step levels (outer and inner) to solve them. 
    111       !! 
    112       !!  The right-hand side (RHS, g1 & g2) are expressed explicitly, and the left-hand side (LHS, f1, f2) implicitly 
    113       !! 
    114       !!  The inner loop makes ~1500 steps max to solve a linear system problem, targetting to solve the implicit terms 
    115       !!  The outer loop makes ~10 steps and updates the explicit terms 
    116       !!  Outer loop uses a modified euler time step, with uC = 0.5 * ( uC + u ) 
     120      !!  The right-hand side (RHS) is explicit 
     121      !!  The left-hand side (LHS) is implicit 
     122      !!  Coriolis is part of explicit terms, whereas ice-ocean drag is implicit 
     123      !! 
     124      !!  Two iteration levels (outer and inner loops) are used to solve the equations 
     125      !! 
     126      !!  The outer loop (OL, typically 10 iterations) is there to deal with the (strong) non-linearities in the equation 
     127      !! 
     128      !!  The inner loop (IL, typically 1500 iterations) is there to solve the linear problem with a line-successive-relaxation algorithm 
     129      !! 
     130      !!  The velocity used in the non-linear terms uses a "modified euler time step" (not sure its the correct term),  
     131      !!! with uk = ( uk-1 + uk-2 ) / 2. 
    117132      !!   
    118133      !!  * Spatial discretization  
     
    143158      !! ** Notes   :  
    144159      !!              
    145       !! References : Zhang and Hibler, JGR 1997; Losch et al., OM 2010.  
     160      !! References : Zhang and Hibler, JGR 1997; Losch et al., OM 2010., Lemieux et al., 2008, 2009, ...  
    146161      !!              
    147162      !!              
    148163      !!------------------------------------------------------------------- 
    149164      !! 
    150 !!!      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step 
    151 !!!      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   ! 
    152 !!!      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
    153 !!!      !! 
    154 !!!      INTEGER ::   ji, jj       ! dummy loop indices 
    155 !!!      INTEGER ::   jter         ! local integers 
    156 !!!      ! 
    157 !!!      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio 
    158 !!!      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling 
    159 !!!      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity 
    160 !!!      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017 
    161 !!!      REAl(wp) ::   zbetau, zbetav 
    162 !!!      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    163 !!!      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars 
    164 !!!      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars 
    165 !!!      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    166 !!!      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast 
    167 !!!      ! 
    168 !!!      REAL(wp) ::   zintb, zintn                                        ! dummy argument 
    169 !!!      REAL(wp) ::   zfac_x, zfac_y 
    170 !!!      REAL(wp) ::   zshear, zdum1, zdum2 
    171 !!!      ! 
    172 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    173 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    174 !!!      ! 
    175 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points 
    176 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points 
    177 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    178 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    179 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
    180 !!!      ! 
    181 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
    182 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
    183 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
    184 !!!      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
    185 !!!      !                                                                 !    ice bottom surface if ice is embedded    
    186 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses 
    187 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points 
    188 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array 
    189 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points 
    190 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points 
     165      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step 
     166      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   ! 
     167      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      ! 
     168      !! 
     169      LOGICAL ::   ll_u_iterate, ll_viterate   ! continue iteration or not 
     170      ! 
     171      INTEGER ::   ji, jj, jn          ! dummy loop indices 
     172      INTEGER ::   jter, i_out, i_inn  !  
     173      INTEGER ::   ji_min, jj_min      ! 
     174      INTEGER ::   nn_zebra_vp         ! number of zebra steps 
     175      INTEGER ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp)       
     176      ! 
     177      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio 
     178      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity 
     179      REAL(wp) ::   zglob_area                                          ! global ice area for diagnostics 
     180      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
     181      REAL(wp) ::   zm2, zm3, zmassU, zmassV                            ! ice/snow mass and volume 
     182      REAL(wp) ::   zdeltat, zdeltat_star, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 
     183      REAL(wp) ::   zp_deltastar_f 
     184      REAL(wp) ::   zfac, zfac1, zfac2, zfac3 
     185      REAL(wp) ::   zt12U, zt11U, zt22U, zt21U, zt122U, zt121U 
     186      REAL(wp) ::   zt12V, zt11V, zt22V, zt21V, zt122V, zt121V 
     187      REAL(wp) ::   zAA3, zw, zuerr_max, zverr_max 
     188      ! 
     189      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice 
     190      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points 
     191      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! Acceleration term contribution to RHS 
     192      REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step 
     193      ! 
     194      REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta at T points 
     195      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zet                        ! Viscosity pre-factors at T points 
     196      REAL(wp), DIMENSION(jpi,jpj) ::   zef                             ! Viscosity pre-factor at F point 
     197      ! 
     198      REAL(wp), DIMENSION(jpi,jpj) ::   zmt                             ! Mass per unit area at t-point 
     199      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! Coriolis factor (m*f) at t-point  
     200      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points 
     201      REAL(wp), DIMENSION(jpi,jpj) ::   zu_c, zv_c                      ! "current" ice velocity (m/s), average of previous two OL iterates 
     202      REAL(wp), DIMENSION(jpi,jpj) ::   zu_cV, zv_cU                    ! "current" u (v) ice velocity at V (U) point (m/s) 
     203      REAL(wp), DIMENSION(jpi,jpj) ::   zu_b, zv_b                      ! velocity at previous sub-iterate 
     204      REAL(wp), DIMENSION(jpi,jpj) ::   zuerr, zverr                    ! absolute U/Vvelocity difference between current and previous sub-iterates 
     205      ! 
     206      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
     207      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
     208      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
     209      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
     210      !                                                                 !    ice bottom surface if ice is embedded    
     211      REAL(wp), DIMENSION(jpi,jpj) ::   zCwU, zCwV                      ! ice-ocean drag pre-factor (rho*c*module(u)) 
     212      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points 
     213      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array 
     214      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points 
     215      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi_rhsu, ztauy_oi_rhsv    ! ice-ocean stress RHS contribution at U-V points 
     216      REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsu, zs2_rhsu, zs12_rhsu   ! internal stress contributions to RHSU 
     217      REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsv, zs2_rhsv, zs12_rhsv   ! internal stress contributions to RHSV 
     218      REAL(wp), DIMENSION(jpi,jpj) ::   zf_rhsu, zf_rhsv                ! U- and V- components of internal force RHS contributions 
     219      REAL(wp), DIMENSION(jpi,jpj) ::   zrhsu, zrhsv                    ! U and V RHS  
     220      REAL(wp), DIMENSION(jpi,jpj) ::   zAU, zBU, zCU, zDU, zEU         ! Linear system coefficients, U equation 
     221      REAL(wp), DIMENSION(jpi,jpj) ::   zAV, zBV, zCV, zDV, zEV, zFV    ! Linear system coefficients, V equation 
     222      REAL(wp), DIMENSION(jpi,jpj) ::   zFU, zFU_prime, zBU_prime       ! Rearranged linear system coefficients, U equation 
     223      REAL(wp), DIMENSION(jpi,jpj) ::   zFV, zFV_prime, zBV_prime       ! Rearranged linear system coefficients, V equation 
    191224!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast) 
    192225!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast) 
    193 !!!      ! 
    194 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays 
    195 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence 
    196  
    197 !!! 
    198 !!!      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
    199 !!!      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small 
    200 !!!      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    201 !!! 
    202 !!!      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    203 !!!      !! --- check convergence 
    204 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice 
    205 !!!      !! --- diags 
    206  
    207          !! --- check convergence 
    208           
    209           
    210          !! --- LSR arrays 
    211          REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                           ! mask at F points for the ice 
    212  
    213          REAL(wp), DIMENSION(jpi,jpj) ::   zaU      , zaV                   ! ice fraction on U/V points 
    214          REAL(wp), DIMENSION(jpi,jpj) ::   zmassU   , zmassV                ! ice mass per unit area at U/V points (kg/m2) 
    215           
    216          REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice_C , zv_ice_C              ! ice velocities averaged between before and mid outer sub-iteration 
    217          REAL(wp), DIMENSION(jpi,jpj) ::   uTmp    , vTmp                   ! uTmp  = velocity at previous *inner* interation  
    218  
    219          REAL(wp), DIMENSION(jpi,jpj) ::   z_etat  , z_zetat                ! VP viscosities at T points 
    220          REAL(wp), DIMENSION(jpi,jpj) ::   z_etaf                           ! VP eta viscosity at F point 
    221           
    222          REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsu, zs2_rhsu, zs12_rhsu    ! Internal stress contribution to RHS of U-equation 
    223          REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsv, zs2_rhsv, zs12_rhsv    ! Internal stress contribution to RHS of V-equation 
    224          REAL(wp), DIMENSION(jpi,jpj) ::   zf_rhsu, zf_rhsv                 ! U- and V- components of internal forces 
    225           
    226 !!!      !! --- SIMIP diags 
     226     ! 
     227      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays 
     228      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence 
     229      ! 
     230      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
     231      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small 
     232      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
     233      !! --- diags 
     234      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3          
     235      !! --- SIMIP diags 
    227236      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
    228237      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) 
     
    231240      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s) 
    232241      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)       
    233                    
    234       !! --- Convergence diags 
    235       REAL(wp), ALLOCATABLE, DIMENSION(:) ::  zdiag_meanke(:) 
    236       REAL(wp), ALLOCATABLE, DIMENSION(:) ::  zdiag_resnor(:) 
    237       REAL(wp), ALLOCATABLE, DIMENSION(:) ::  zdiag_veldif(:) 
    238        
     242                         
    239243      !!---------------------------------------------------------------------------------------------------------------------- 
    240244 
    241245      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' 
    242        
    243       ! Namelist parameters to outsource (once we all are super happy) 
    244       nn_nouter_vp      = 3            ! Number of outer iterations (mitGCM 10)  
    245       nn_ninner_vp      = 1500         ! Number of inner-loop iterations (mitGCM 1500) 
    246       rn_deltamin       = 2.0e-9       ! minimum delta value following Lemieux and Dupont, GMD 2020                           
    247       ln_visc_VP_smooth = .FALSE.      ! zeta viscosity smoothing (if yes, tanh function of Lemieux and Tremblay JGR 2009) 
    248  
    249       ! Linear system parameters 
    250       ln_zebra          = .FALSE.      ! do two passages for the linear system problem, one every odd j-band, one every even one (recommended by Martin Losch) 
    251       IF ( ln_zebra ) THEN; nn_zebra_step = 2; ELSE; nn_zebra_step = 1; ENDIF  
    252  
    253       rn_relaxfacU = 0.95 ! U-relaxation factor for linear problem (Losch 0.95, Zhang 1.8) 
    254       rn_relaxfacV = 0.95 ! V-relaxation factor    
    255        
    256       rn_uerr_max = 0.8_wp ! maximum error on velocity, above which a forcing error is considered 
    257       rn_uerr_min = 1.0d-4 ! velocity difference beyond which the linear system procedure is stopped 
    258       nn_lsr_cvgcheck = 5  ! iterations every each convergence is checked 
    259  
     246             
    260247      !------------------------------------------------------------------------------! 
    261       ! 1) Initialization 
     248      ! 
     249      ! --- Initialization 
     250      ! 
    262251      !------------------------------------------------------------------------------! 
     252       
     253      ! for diagnostics and convergence tests 
     254      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 
     255      DO jj = 1, jpj 
     256         DO ji = 1, jpi 
     257            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice 
     258            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 
     259         END DO 
     260      END DO 
     261       
     262      IF ( ln_zebra_vp ) THEN; nn_nzebra_vp = 2; ELSE; nn_nzebra_vp = 1; ENDIF  
     263       
     264      nn_nvp = nn_out_vp * nn_inn_vp ! maximum number of iterations 
     265       
    263266      zrhoco = rau0 * rn_cio  
    264267 
     
    272275      zs12(:,:) = pstress12_i(:,:) 
    273276       
    274       ! Convergence checks 
    275       IF( ln_rhg_chkcvg ) THEN  
    276        
     277      ! Initialise convergence checks 
     278      IF( ln_rhg_chkcvg ) THEN 
     279       
     280         ! ice area for global mean kinetic energy 
    277281         zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) ) ! global ice area (km2) 
    278          z1_32      = 1._wp / 32._wp 
    279  
    280          ! --- mean kinetic energy       
    281          ALLOCATE( zdiag_meanke(nn_outer_vp*nn_inner_vp) )  ! mean kinetic energy           
    282          ALLOCATE( zdiag_resnor(nn_outer_vp*nn_inner_vp) )  ! residual norm of the linear system   
    283          ALLOCATE( zdiag_u_idif(nn_outer_vp*nn_inner_vp) )  ! velocity difference between two iterations 
    284          ALLOCATE( zdiag_v_idif(nn_outer_vp*nn_inner_vp) )  ! velocity difference between two iterations 
    285           
    286          zdiag_resnor(:) = 0._wp 
    287          zdiag_meanke(:) = 0._wp 
    288          zdiag_u_idif(:) = 0._wp 
    289          zdiag_v_idif(:) = 0._wp 
    290282          
    291283      ENDIF 
    292284 
    293       ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
     285      ! Landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
     286      ! MV: Not working yet... 
    294287      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile 
    295288      ELSE                         ;   zkt = 0._wp 
     
    297290 
    298291      !------------------------------------------------------------------------------! 
    299       ! 2) Time-independent quantities 
     292      ! 
     293      ! --- Time-independent quantities 
     294      ! 
    300295      !------------------------------------------------------------------------------! 
    301296       
    302297      CALL ice_strength ! strength at T points 
    303298       
    304       !-------- 
    305       ! F-mask       
    306       !-------- 
     299      !------------------------------ 
     300      ! -- F-mask       (code from EVP) 
     301      !------------------------------ 
    307302      ! MartinV:  
    308303      ! In EVP routine, zfmask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip) 
    309       ! I don't know exactly where it should be applied in here 
    310        
    311       ! ocean/land mask 
    312       DO jj = 1, jpjm1 
    313          DO ji = 1, jpim1      ! NO vector opt. 
     304      ! I am not sure the same recipe applies here 
     305       
     306      ! - ocean/land mask 
     307      DO jj = 1, jpj - 1 
     308         DO ji = 1, jpi - 1 
    314309            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 
    315310         END DO 
     
    318313 
    319314      ! Lateral boundary conditions on velocity (modify zfmask) 
    320       DO jj = 2, jpjm1 
    321          DO ji = fs_2, fs_jpim1   ! vector opt. 
     315      ! Can be computed once for all, at first time step, for all rheologies 
     316      DO jj = 2, jpj - 1 
     317         DO ji = 2, jpi - 1   ! vector opt. 
    322318            IF( zfmask(ji,jj) == 0._wp ) THEN 
    323319               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 
     
    326322         END DO 
    327323      END DO 
    328       DO jj = 2, jpjm1 
     324      DO jj = 2, jpj - 1 
    329325         IF( zfmask(1,jj) == 0._wp ) THEN 
    330326            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 
    331327         ENDIF 
    332328         IF( zfmask(jpi,jj) == 0._wp ) THEN 
    333             zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 
     329            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpi - 1,jj,1), umask(jpi,jj-1,1) ) ) 
    334330        ENDIF 
    335331      END DO 
    336       DO ji = 2, jpim1 
     332      DO ji = 2, jpi - 1 
    337333         IF( zfmask(ji,1) == 0._wp ) THEN 
    338334            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 
    339335         ENDIF 
    340336         IF( zfmask(ji,jpj) == 0._wp ) THEN 
    341             zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 
     337            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpj - 1,1) ) ) 
    342338         ENDIF 
    343339      END DO 
     
    345341      CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp ) 
    346342       
    347       !--------------------------------------------------------------------------------------------------------- 
    348       ! 3) Time-independent contributions of acceleration, ocean drag, coriolis, atmospheric drag, surface tilt 
    349       !--------------------------------------------------------------------------------------------------------- 
    350       !  
    351       ! --- Compute all terms & factors independent of velocities, or only depending on velocities at previous time step 
    352        
    353       DO jj = 2, jpjm1 
    354          DO ji = fs_2, fs_jpim1 
     343      !---------------------------------------------------------------------------------------------------------- 
     344      ! -- Time-independent pre-factors for acceleration, ocean drag, coriolis, atmospheric drag, surface tilt 
     345      !---------------------------------------------------------------------------------------------------------- 
     346      ! Compute all terms & factors independent of velocities, or only depending on velocities at previous time step       
     347       
     348      ! sea surface height 
     349      !    embedded sea ice: compute representative ice top surface 
     350      !    non-embedded sea ice: use ocean surface for slope calculation 
     351      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
     352       
     353      DO jj = 2, jpj - 1 
     354         DO ji = 2, jpi - 1 
    355355 
    356356            ! Ice fraction at U-V points 
     
    359359 
    360360            ! Snow and ice mass at U-V points 
    361             zm1             = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )  
     361            zmt(ji,jj)      = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )  
    362362            zm2             = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) ) 
    363363            zm3             = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) ) 
    364364             
    365             zmassU          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    366             zmassV          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     365            zmassU          = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
     366            zmassV          = 0.5_wp * ( zmt(ji,jj) * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    367367                           
    368             ! Mass per unit time 
     368            ! Mass per unit area divided by time step 
    369369            zmassU_t(ji,jj) = zmassU * r1_rdtice 
    370370            zmassV_t(ji,jj) = zmassV * r1_rdtice 
     
    379379 
    380380            ! Coriolis factor at T points (m*f) 
    381             zmf(ji,jj)      = zm1 * ff_t(ji,jj) 
     381            zmf(ji,jj)      = zmt(ji,jj) * ff_t(ji,jj) 
    382382             
    383383            ! Wind stress 
     
    393393            zmsk00y(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
    394394 
    395             ! switches 
     395            ! switches  
    396396            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
    397397            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
     
    403403             
    404404      !------------------------------------------------------------------------------! 
    405       ! 4) Start outer loop, update velocities 
     405      ! 
     406      ! --- Start outer loop 
     407      ! 
    406408      !------------------------------------------------------------------------------! 
    407409       
    408       zu_ice_C(:,:) = u_ice(:,:) 
    409       zv_ice_C(:,:) = v_ice(:,:) 
    410        
    411       i_iter = 0 
    412  
    413       DO i_out = 1, nn_nouter_vp 
    414        
    415          i_iter = i_iter + 1 
    416        
    417          ! EVP CODE 
    418          ! l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
    419          ! 
    420          ! convergence test 
    421          IF( ln_rhg_chkcvg ) THEN 
    422             DO jj = 1, jpj 
    423                DO ji = 1, jpi 
    424                   zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 
    425                   zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 
    426                END DO 
    427             END DO 
    428          ENDIF 
    429           
    430          ! Zhang and Hibler 97 time stepping 
    431          ! identified as the most rapid test with mitGCM and recommendation by Martin Losch 
    432          zu_ice_C(:,:) = 0.5_wp * ( u_ice(:,:) + zu_ice_C(:,:) ) 
    433          zv_ice_C(:,:) = 0.5_wp * ( v_ice(:,:) + zv_ice_C(:,:) ) 
    434           
    435       !------------------------------------------------------------------------------! 
    436       ! 5) Right-hand side of linear problem 
    437       !------------------------------------------------------------------------------! 
    438  
    439          !       
     410      zu_c(:,:) = u_ice(:,:) 
     411      zv_c(:,:) = v_ice(:,:) 
     412       
     413      jter = 0 
     414 
     415      DO i_out = 1, nn_nout_vp 
     416                
     417         ! 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}) 
     419         ! Also used in Hibler and Ackley (1983); Zhang and Hibler (1997); Lemieux and Tremblay (2009) 
     420         zu_c(:,:) = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) ) 
     421         zv_c(:,:) = 0.5_wp * ( v_ice(:,:) + zv_c(:,:) ) 
     422          
     423         !------------------------------------------------------------------------------! 
     424         ! 
     425         ! --- Right-hand side (RHS) of the linear problem 
     426         ! 
     427         !------------------------------------------------------------------------------! 
    440428         ! In the outer loop, one needs to update all RHS terms 
    441429         ! with explicit velocity dependencies (viscosities, coriolis, ocean stress) 
    442          ! 
    443        
    444          !-------------------------------------------------------- 
    445          ! 5.1 Strain rates, viscosities and replacement pressure 
    446          !-------------------------------------------------------- 
     430         ! as a function of uc 
     431         ! 
     432       
     433         !------------------------------------------ 
     434         ! -- Strain rates, viscosities and P/Delta 
     435         !------------------------------------------ 
    447436 
    448437         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
    449          DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
    450             DO ji = 1, jpim1 
     438         DO jj = 1, jpj - 1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
     439            DO ji = 1, jpi - 1 
    451440 
    452441               ! shear at F points 
    453                zds(ji,jj) = ( ( zu_ice_C(ji,jj+1) * r1_e1u(ji,jj+1) - zu_ice_C(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    454                   &         + ( zv_ice_C(ji+1,jj) * r1_e2v(ji+1,jj) - zv_ice_C(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     442               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)   & 
     443                  &         + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    455444                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 
    456445 
     
    458447         END DO 
    459448          
    460          CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) 
    461  
    462          DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
    463             DO ji = 2, jpi ! no vector loop 
     449         CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan 
     450 
     451         DO jj = 2, jpj - 1    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
     452            DO ji = 2, jpi - 1 !  
    464453 
    465454               ! shear**2 at T points (doc eq. A16) 
     
    469458               
    470459               ! divergence at T points 
    471                zdiv  = ( e2u(ji,jj) * zu_ice_C(ji,jj) - e2u(ji-1,jj) * zu_ice_C(ji-1,jj)   & 
    472                   &    + e1v(ji,jj) * zv_ice_C(ji,jj) - e1v(ji,jj-1) * zv_ice_C(ji,jj-1)   & 
     460               zdiv  = ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj)   & 
     461                  &    + e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1)   & 
    473462                  &    ) * r1_e1e2t(ji,jj) 
    474463               zdiv2 = zdiv * zdiv 
    475464                
    476465               ! tension at T points 
    477                zdt  = ( ( zu_ice_C(ji,jj) * r1_e2u(ji,jj) - zu_ice_C(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    478                   &   - ( zv_ice_C(ji,jj) * r1_e1v(ji,jj) - zv_ice_C(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     466               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)   & 
     467                  &   - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    479468                  &   ) * r1_e1e2t(ji,jj) 
    480469               zdt2 = zdt * zdt 
     
    486475               zdeltat_star = MAX( zdeltat, rn_delta_min ) 
    487476                
    488 !               IF ( ln_visc_VP_smooth ) THEN ! to code, other regularization from Lemieux and Tremblay (2009) 
    489 !                   zdeltat_star = ... 
    490 !               ENDIF 
     477               IF ( ln_smooth_vp ) zdelta_star = zdeltat + rn_delta_min 
    491478                
    492479               ! P/delta at T-points 
     
    494481                
    495482               ! Temporary zzt and zet factors at T-points 
    496                zzt(ji,jj)     = zp_deltastar_t * r1_e1e2t(ji,jj) 
     483               zzt(ji,jj)     = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj) 
    497484               zet(ji,jj)     = zzt(ji,jj)     * z1_ecc2  
    498485                           
     
    500487         END DO 
    501488 
    502          CALL lbc_lnk( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. ) 
    503          CALL lbc_lnk( 'icedyn_rhg_vp', zzt , 'T', 1. ) 
    504          CALL lbc_lnk( 'icedyn_rhg_vp', zet, 'T', 1. ) 
    505  
    506          DO jj = 1, jpjm1 
    507             DO ji = 1, jpim1 
     489         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 
     490 
     491         DO jj = 1, jpj - 1 
     492            DO ji = 1, jpi - 1 
    508493                
    509494               ! P/delta* at F points 
     
    519504          
    520505         !--------------------------------------------------- 
    521          ! 5.2 Ocean-ice drag and Coriolis RHS contributions 
     506         ! -- Ocean-ice drag and Coriolis RHS contributions 
    522507         !--------------------------------------------------- 
    523508          
    524          DO jj = 2, jpjm1 
    525              DO ji = fs_2, fs_jpim1 
     509         DO jj = 2, jpj - 1 
     510             DO ji = 2, jpi - 1 
    526511                 
    527512                !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation) 
    528                 zu_ice_CV(ji,jj)     = 0.25_wp * ( zu_ice_C(ji,jj) + zu_ice_C(ji-1,jj) + zu_ice_C(ji,jj+1) + zu_ice_C(ji-1,jj+1) ) * vmask(ji,jj,1) 
    529                 zv_ice_CU(ji,jj)     = 0.25_wp * ( zv_ice_C(ji,jj) + zv_ice_C(ji,jj-1) + zv_ice_C(ji+1,jj) + zv_ice_C(ji+1,jj-1) ) * umask(ji,jj,1) 
    530  
     513                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) 
     514                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) 
     515                 
    531516                !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3) 
    532                 zCwU(ji,jj)          = zaU(ji,jj) * zrhoco * SQRT( ( zu_ice_C (ji,jj) - u_oce (ji,jj) ) * ( zu_ice_C (ji,jj) - u_oce (ji,jj) )  & 
    533                   &                                              + ( zv_ice_CU(ji,jj) - v_oceU(ji,jj) ) * ( zv_ice_CU(ji,jj) - v_oceU(ji,jj) ) ) 
    534                 zCwV(ji,jj)          = zaV(ji,jj) * zrhoco * SQRT( ( zv_ice_C (ji,jj) - v_oce (ji,jj) ) * ( zv_ice_C (ji,jj) - v_oce (ji,jj) )  & 
    535                   &                                              + ( zu_ice_CV(ji,jj) - u_oceV(ji,jj) ) * ( zu_ice_CV(ji,jj) - u_oceV(ji,jj) ) ) 
     517                zCwU(ji,jj)          = zaU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) )  & 
     518                  &                                              + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) ) 
     519                zCwV(ji,jj)          = zaV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) )  & 
     520                  &                                              + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) ) 
    536521                  
    537522                !--- Ocean-ice drag contributions to RHS  
     
    542527                ! Note Lemieux et al 2008 recommend to do that implicitly, but I don't really see how this could be done 
    543528                zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  & 
    544                            &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * zv_ice_C(ji  ,jj) + e1v(ji  ,jj-1) * zv_ice_C(ji  ,jj-1) )  & 
    545                            &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_ice_C(ji+1,jj) + e1v(ji+1,jj-1) * zv_ice_C(ji+1,jj-1) ) ) 
     529                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * zv_c(ji  ,jj) + e1v(ji  ,jj-1) * zv_c(ji  ,jj-1) )  & 
     530                           &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_c(ji+1,jj) + e1v(ji+1,jj-1) * zv_c(ji+1,jj-1) ) ) 
    546531                            
    547532                ! --- V-component of Coriolis Force (energy conserving formulation) 
    548533                zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  & 
    549                            &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * zu_ice_C(ji,jj  ) + e2u(ji-1,jj  ) * zu_ice_C(ji-1,jj  ) )  & 
    550                            &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_ice_C(ji,jj+1) + e2u(ji-1,jj+1) * zu_ice_C(ji-1,jj+1) ) ) 
     534                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * zu_c(ji,jj  ) + e2u(ji-1,jj  ) * zu_c(ji-1,jj  ) )  & 
     535                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) ) 
    551536          
    552537             END DO 
    553538         END DO 
    554539          
    555          ! ??? lbclnk ??? 
    556  
    557          !-------------------------------------- 
    558          ! 5.3 Internal stress RHS contribution 
    559          !-------------------------------------- 
    560  
    561          ! --- Stress contributions at t-points          
     540         ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system,  
     541         ! so there is no need for lbclnk on drag and coriolis 
     542 
     543         !------------------------------------- 
     544         ! -- Internal stress RHS contribution 
     545         !------------------------------------- 
     546 
     547         ! --- Stress contributions at T-points          
    562548         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
    563             DO ji = 2, jpi ! no vector loop 
     549            DO ji = 2, jpi !  
    564550             
    565551               ! sig1 contribution to RHS of U-equation at T-points 
    566                zs1_rhsu(ji,jj) =   zzt(ji,jj) * ( e1v(ji,jj)    * zv_ice_C(ji,jj) - e1v(ji,jj-1)    * zv_ice_C(ji,jj-1) - 1.0_wp ) 
     552               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 ) 
    567553                                             
    568554               ! sig2 contribution to RHS of U-equation at T-points             
    569                zs2_rhsu(ji,jj) = - zet(ji,jj) * ( r1_e1v(ji,jj) * zv_ice_C(ji,jj) - r1_e1v(ji,jj-1) * zv_ice_C(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  
     555               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)  
    570556 
    571557               ! sig1 contribution to RHS of V-equation at T-points 
    572                zs1_rhsv(ji,jj) =   zzt(ji,jj) * ( e2u(ji,jj) * zu_ice_C(ji,jj)    - e2u(ji-1,jj)    * zu_ice_C(ji-1,jj) - 1.0_wp ) 
     558               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 ) 
    573559 
    574560               ! sig2 contribution to RHS of V-equation  at T-points 
    575                zs2_rhsv(ji,jj) =   zet(ji,jj) * ( r1_e2u(ji,jj) * zu_ice_C(ji,jj) - r1_e2u(ji-1,jj) * zu_ice_C(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) 
     561               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) 
    576562 
    577563            END DO 
    578564         END DO 
    579565          
    580          ! ??? lbclnk ??? 
     566         ! a priori, no lbclnk, because rhsu is only used in the inner domain 
    581567          
    582568         ! --- Stress contributions at f-points          
    583          ! Note from MartinV: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing 
     569         ! MV NOTE: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing 
    584570         ! My guess is that this is the way to enforce boundary conditions on strain rate tensor 
    585571 
    586          DO jj = 1, jpjm1 
    587             DO ji = 1, jpim1 
     572         DO jj = 1, jpj - 1 
     573            DO ji = 1, jpi - 1 
    588574                
    589575               ! sig12 contribution to RHS of U equation at F-points  
    590                zs12_rhsu(ji,jj) = - zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_ice_C(ji+1,jj) - r1_e2v(ji,jj) * zv_ice_C(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj) 
     576               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) * zfmask(ji,jj) 
    591577                
    592578               ! sig12 contribution to RHS of V equation at F-points 
    593                zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_ice_C(ji,jj+1) - r1_e1u(ji,jj) * zu_ice_C(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj) 
     579               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) * zfmask(ji,jj) 
    594580 
    595581            END DO 
    596582         END DO 
    597583          
    598          ! ??? lbcblnk 
    599  
    600          ! --- Internal forces contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002) 
    601          DO jj = 2, jpjm1 
    602             DO ji = fs_2, fs_jpim1                
     584         ! a priori, no lbclnk, because rhsu are only used in the inner domain 
     585 
     586         ! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002) 
     587         ! OPT: merge with next loop and use intermediate scalars for zf_rhsu 
     588          
     589         DO jj = 2, jpj - 1 
     590            DO ji = 2, jpi - 1                
    603591               ! --- U component of internal force contribution to RHS at U points 
    604592               zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * &  
     
    606594                  &      +         r1_e2u(ji,jj) * ( e2t(ji+1,jj) * e2t(ji+1,jj) * zs2_rhsu(ji+1,jj) - e2t(ji,jj)   * e2t(ji,jj)   * zs2_rhsu(ji,jj) )     & 
    607595                  &      + 2._wp * r1_e1u(ji,jj) * ( e1f(ji,jj)   * e1f(ji,jj)   * zs12_rhsu(ji,jj)  - e1f(ji,jj-1) * e1f(ji,jj-1) * zs12_rhsu(ji,jj-1) )  
     596                   
    608597               ! --- V component of internal force contribution to RHS at V points 
    609598               zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * & 
     
    611600                  &      +         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) )     & 
    612601                  &      + 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) ) 
     602                   
    613603            END DO 
    614604         END DO 
    615605          
    616          ! ??? lbclnk ??? 
    617  
    618606         !--------------------------- 
    619          ! 5.4 Sum RHS contributions 
     607         ! -- Sum RHS contributions 
    620608         !--------------------------- 
    621          DO jj = 2, jpjm1 
    622             DO ji = fs_2, fs_jpim1 
     609         ! 
     610         ! OPT: could use intermediate scalars to reduce memory access 
     611         DO jj = 2, jpj - 1 
     612            DO ji = 2, jpi - 1 
    623613             
    624614               ! still miss ice ocean stress and acceleration contribution 
     
    629619         END DO 
    630620          
    631       ! ??? lbclnk ??? 
     621         ! inner domain calculations -> no lbclnk 
    632622      
    633       !---------------------------------------------------------------------------------------! 
    634       ! 6) Linear system matrix 
    635       !---------------------------------------------------------------------------------------! 
     623         !---------------------------------------------------------------------------------------! 
     624         ! 
     625         ! --- Linear system matrix 
     626         ! 
     627         !---------------------------------------------------------------------------------------! 
    636628       
    637629         ! Linear system matrix contains all implicit contributions  
     
    641633         ! AU * u_{i-1,j} + BU * u_{i,j}   + CU * u_{i+1,j} 
    642634         !                = DU * u_{i,j-1} + EU * u_{i,j+1} + RHS          (! my convention, not the same as ZH97 )          
    643          ! 
    644          !-------------------------------------- 
    645          ! 6.1 Internal forces LHS contribution 
    646          !-------------------------------------- 
    647          ! 
    648          ! Note 1: martin losch applies boundary condition to BU in mitGCM - check whether it is necessary here ? 
    649          ! Note 2: "T" factor calculations can be optimized by putting things out of the loop  
     635          
     636         ! MV Note 1: martin losch applies boundary condition to BU in mitGCM - check whether it is necessary here ? 
     637         ! MV Note 2: "T" factor calculations can be optimized by putting things out of the loop  
    650638         !         only zzt and zet are iteration-dependent, other only depend on scale factors 
    651          ! 
    652          ! --- U-component 
    653          ! 
    654          ! "T" factors (intermediate results) 
    655          ! 
    656          zfac       = 0.5_wp * r1_e1e2u(ji,jj) 
    657          zfac1      =         zfac * e2u(ji,jj) 
    658          zfac2      =         zfac * r1_e2u(ji,jj) 
    659          zfac3      = 2._wp * zfac * r1_e1u(ji,jj) 
    660  
    661          zt12U      = - zfac1 * zzt(ji+1,jj) 
    662          zt11U      =   zfac1 * zzt(ji,jj) 
    663           
    664          zt22U      = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 
    665          zt21U      =   zfac2 * zet(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj) 
    666           
    667          zt122U     = - zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj) 
    668          zt121U     =   zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 
    669  
    670          ! 
    671          ! Linear system coefficients 
    672          ! 
    673          zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj) 
    674          zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj) 
    675          zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj) 
    676           
    677          zDU(ji,jj) =   zt121U * r1_e1u(ji,jj-1) 
    678          zEU(ji,jj) =   zt122U * r1_e1u(ji,jj+1) 
    679           
    680          ! 
    681          ! --- V-component 
    682          ! 
    683          ! "T" factors (intermediate results) 
    684          ! 
    685          zfac       = 0.5_wp * r1_e1e2v(ji,jj) 
    686          zfac1      =         zfac * e2v(ji,jj) 
    687          zfac2      =         zfac * r1_e1v(ji,jj) 
    688          zfac3      = 2._wp * zfac * r1_e2v(ji,jj) 
    689           
    690          zt12V      = - zfac1 * zzt(ji,jj+1) 
    691          zt11V      =   zfac1 * zzt(ji,jj) 
    692           
    693          zt22V      =   zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) 
    694          zt21V      = - zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj) 
    695           
    696          zt122V     =   zfac3 * zef(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj) 
    697          zt121V     = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 
    698           
    699          ! 
    700          ! Linear system coefficients 
    701          ! 
    702          zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1) 
    703          zBV(ji,jj) =  ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj) 
    704          zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1) 
    705           
    706          zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here 
    707          zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj) 
    708                    
    709          !--------------------------------- 
    710          ! Ocean-ice drag LHS contribution 
    711          !--------------------------------- 
    712          zBU(i,j) = zBU(i,j) + zCwU(ji,jj) 
    713  
    714          !------------------------------- 
    715          ! Acceleration LHS contribution 
    716          !------------------------------- 
    717          zBU(i,j) = zBU(i,j) + zmassU_t(ji,jj)  
     639                   
     640         DO ji = 2, jpi - 1 ! internal domain do loop 
     641            DO jj = 2, jpj - 1 
     642 
     643               !------------------------------------- 
     644               ! -- Internal forces LHS contribution 
     645               !------------------------------------- 
     646               ! 
     647               ! --- U-component 
     648               ! 
     649               ! "T" factors (intermediate results) 
     650               ! 
     651               zfac       = 0.5_wp * r1_e1e2u(ji,jj) 
     652               zfac1      =         zfac * e2u(ji,jj) 
     653               zfac2      =         zfac * r1_e2u(ji,jj) 
     654               zfac3      = 2._wp * zfac * r1_e1u(ji,jj) 
     655 
     656               zt12U      = - zfac1 * zzt(ji+1,jj) 
     657               zt11U      =   zfac1 * zzt(ji,jj) 
     658          
     659               zt22U      = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 
     660               zt21U      =   zfac2 * zet(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj) 
     661          
     662               zt122U     = - zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj) 
     663               zt121U     =   zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 
     664                
     665               ! 
     666               ! Linear system coefficients 
     667               ! 
     668               zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj) 
     669               zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj) 
     670               zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj) 
     671          
     672               zDU(ji,jj) =   zt121U * r1_e1u(ji,jj-1) 
     673               zEU(ji,jj) =   zt122U * r1_e1u(ji,jj+1) 
     674               
     675               ! 
     676               ! --- V-component 
     677               ! 
     678               ! "T" factors (intermediate results) 
     679               ! 
     680               zfac       = 0.5_wp * r1_e1e2v(ji,jj) 
     681               zfac1      =         zfac * e2v(ji,jj) 
     682               zfac2      =         zfac * r1_e1v(ji,jj) 
     683               zfac3      = 2._wp * zfac * r1_e2v(ji,jj) 
     684          
     685               zt12V      = - zfac1 * zzt(ji,jj+1) 
     686               zt11V      =   zfac1 * zzt(ji,jj) 
     687          
     688               zt22V      =   zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) 
     689               zt21V      = - zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj) 
     690          
     691               zt122V     =   zfac3 * zef(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj) 
     692               zt121V     = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 
     693          
     694               ! 
     695               ! Linear system coefficients 
     696               ! 
     697               zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1) 
     698               zBV(ji,jj) =  ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj) 
     699               zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1) 
     700          
     701               zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here 
     702               zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj) 
     703                   
     704               !----------------------------------------------------- 
     705               ! -- Ocean-ice drag and acceleration LHS contribution 
     706               !----------------------------------------------------- 
     707               zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj) 
     708               zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 
     709          
     710            END DO 
     711         END DO 
    718712                
    719713      !------------------------------------------------------------------------------! 
    720       ! 6) inner loop: solve linear system, check convergence 
     714      ! 
     715      ! --- Inner loop: solve linear system, check convergence 
     716      ! 
    721717      !------------------------------------------------------------------------------! 
    722718                
    723719         ! Inner loop solves the linear problem .. requires 1500 iterations 
    724          doIterate4u = .TRUE. 
    725          doIterate4v = .TRUE. 
    726  
    727          DO i_inner = 1, nn_ninner_vp ! inner loop iterations 
    728           
    729             ! initialize residual ... 
    730              
    731             ... 
    732  
    733             IF ( doIterate4u .OR. doIterate4v ) THEN 
    734  
    735                zu_p(:,:) = u_ice(:,:) ! previous sub-iteration value 
    736                zv_p(:,:) = v_ice(:,:) 
    737  
    738                                           ! ---------------------------- ! 
    739                IF ( doIterate4u ) THEN    ! --- Solve for u-velocity --- ! 
    740                                           ! ---------------------------- ! 
     720         ll_u_iterate = .TRUE. 
     721         ll_v_iterate = .TRUE. 
     722 
     723         DO i_inn = 1, nn_ninn_vp ! inner loop iterations 
     724          
     725            ! mitgcm computes initial value of residual 
     726            jter = jter + 1 
     727            l_full_nf_update = jter == nn_nvp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
     728 
     729            IF ( ll_u_iterate .OR. ll_v_iterate ) THEN 
     730             
     731 
     732               zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 
     733               zv_b(:,:) = v_ice(:,:) 
     734 
     735                                           ! ---------------------------- ! 
     736               IF ( ll_u_iterate ) THEN    ! --- Solve for u-velocity --- ! 
     737                                           ! ---------------------------- ! 
    741738       
    742739                  ! what follows could be subroutinized... 
    743740                   
    744                   DO k = 1, nn_zebra_step ! "zebra" loop 
    745        
    746                      IF ( k == 1 ); THEN; j_min = 2; ELSE; j_min = 3 
    747           
    748                      zFU(:,:) = 0._wp 
    749           
    750                      DO jj = j_min, jpjm1, nn_zebra_step 
     741                  DO jn = 1, nn_nzebra_vp ! "zebra" loop (! red-black-sor!!! ) 
     742                   
     743                     ! OPT: could be even better optimized with a true red-black SOR 
     744       
     745                     IF ( jn == 1 ) THEN   ;   jj_min = 2  
     746                     ELSE                  ;   jj_min = 3 
     747                     ENDIF 
     748          
     749                     zFU(:,:)       = 0._wp 
     750                     zFU_prime(:,:) = 0._wp 
     751                     zBU_prime(:,:) = 0._wp 
     752          
     753                     DO jj = jj_min, jpj - 1, nn_nzebra_vp 
    751754 
    752755                        !------------------------------------------- 
    753                         ! 1) Tridiagonal system solver for each row 
     756                        ! -- Tridiagonal system solver for each row 
    754757                        !------------------------------------------- 
     758                        ! 
     759                        ! MV - I am in doubts whether the way I coded it is reproducible - ask Gurvan 
     760                        ! 
    755761                        ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F 
    756762 
     763                        ! - Right-hand side of tridiagonal system (zFU) 
     764                        DO ji = 2, jpi - 1     
     765 
     766                           ! boundary condition substitution 
     767                           ! see Zhang and Hibler, 1997, Appendix B 
     768                           ! MV NOTE possibly not fully appropriate 
     769                           zAA3 = 0._wp 
     770                           IF ( ji == 2 )          zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj) 
     771                           IF ( ji == jpi - 1 )    zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj) 
     772      
     773                           ! right hand side 
     774                           zFU(ji,jj) = ( zrhsu(ji,jj) &                                 ! right-hand side terms 
     775                               &      +   zAA3                                           ! boundary condition translation 
     776                               &      +   zDU(ji,jj) * u_ice(ji,jj-1)                    ! internal force, j-1 
     777                               &      +   zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1 
     778                     
     779                        END DO 
     780                         
     781                        ! - Thomas Algorithm  
     782                        ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 
     783                        ! Forward sweep 
     784                        DO ji = 3, jpi - 1 
     785                           zw         = zAU(ji,jj) / zBU(ji-1,jj) 
     786                           zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 
     787                           zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 
     788                        END DO 
     789             
     790                        ! Backward sweep 
     791                        ! MV I don't see how this will be reproducible 
     792                        u_ice(jpi-1,jj)     = zFU_prime(jpi-1,jj) / zBU_prime(jpi-1,jj) * umask(jpi-1,jj,1)                     ! do last row first 
     793                        DO ji = jpi-2 , 2, -1 ! all other rows    !   MV OPT: could be turned into forward loop (by substituting ji) 
     794                           u_ice(ji,jj)    = zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji,jj+1) * umask(ji,jj,1) / zBU_prime(ji,jj) ! 
     795                        END DO             
     796             
     797                        !---------------        
     798                        ! -- Relaxation 
     799                        !--------------- 
     800                        DO ji = 2, jpi - 1     
     801                           u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) 
     802                        END DO 
     803 
     804                     END DO ! jj 
     805 
     806                  END DO ! zebra loop 
     807                   
     808               ENDIF !   ll_u_iterate 
     809                
     810               !                           ! ---------------------------- ! 
     811               IF ( ll_v_iterate ) THEN    ! --- Solve for V-velocity --- ! 
     812            !                              ! ---------------------------- ! 
     813                                           
     814                  ! MV OPT: what follows could be subroutinized... 
     815                   
     816                  DO jn = 1, nn_nzebra_vp ! "zebra" loop 
     817       
     818                     IF ( jn == 1 ) THEN   ;   ji_min = 2  
     819                     ELSE                  ;   ji_min = 3 
     820                     ENDIF 
     821          
     822                     zFV(:,:)       = 0._wp 
     823                     zFV_prime(:,:) = 0._wp 
     824                     zBV_prime(:,:) = 0._wp 
     825 
     826                     DO ji = ji_min, jpi - 1, nn_nzebra_vp  
     827                         
     828                        !!! It is intentional to have a ji then jj loop for V-velocity 
     829                        !!! ZH97 explain it is critical for convergence speed 
     830 
     831                        !------------------------------------------- 
     832                        ! -- Tridiagonal system solver for each row 
     833                        !------------------------------------------- 
     834                        ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F 
     835 
    757836                        ! --- Right-hand side of tridiagonal system (zFU) 
    758                         DO ji = fs_2, fs_jpim1     
    759  
    760                            ! boundary condition substitution (check whether is correctly applied !!!) 
     837                        DO jj = 2, jpj - 1 
     838 
     839                           ! boundary condition substitution (check it is correctly applied !!!) 
    761840                           ! see Zhang and Hibler, 1997, Appendix B 
    762841                           zAA3 = 0._wp 
    763                            IF ( ji .EQ. fs_2 )     zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj) 
    764                            IF ( ji .EQ. fs_jpim1 ) zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj) 
     842                           IF ( jj .EQ. 2 )       zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1) 
     843                           IF ( jj .EQ. jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1) 
    765844      
    766845                           ! right hand side 
    767                            zFU(ji,jj) = zrhsu(ji,jj) &               ! right-hand side terms 
    768                                &      + zAA3                        ! boundary condition translation 
    769                                &      + zDU(ji,jj) * u_ice(ji,jj-1) ! internal force, j-1 
    770                                &      + zEU(ji,jj) * u_ice(ji,jj+1) ! internal force, j+1 
     846                           zFV(ji,jj) = ( zrhsv(ji,jj) &              ! right-hand side terms 
     847                               &        + zAA3                        ! boundary condition translation 
     848                               &        + zDV(ji,jj) * v_ice(ji-1,jj) ! internal force, j-1 
     849                               &        + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1) ! internal force, j+1 
    771850                     
    772                            zFU(ji,jj) = zFU(ji,jj) * umask(ji,jj)     ! mask - useful ? 
    773          
    774851                        END DO 
    775852        
     
    777854                        ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 
    778855                        ! Forward sweep 
    779                         DO ji = fs_2 + 1, fs_jpim1 
    780                            zw         = zAU(ji,jj) / zBU(ji-1,jj) 
    781                            zBU(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 
    782                            zFU(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 
     856                        DO jj = 3, jpj - 1 
     857                           zw         = zAV(ji,jj) / zBV(ji,jj-1) 
     858                           zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 
     859                           zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1)  
    783860                        END DO 
    784861             
    785862                        ! Backward sweep 
    786                         u_ice(fs_jpim1,jj) = zFU(ji) / zBU(ji) ! last row 
    787                         DO ji  = fs_jpim - 1, fs_2, -1 ! can be turned into forward row by substituting ji if needed 
    788                            u_ice(ji,jj)    = zFU(ji) - zCU(ji,jj) * u_ice(ji,jj+1) / zBU(ji,jj) 
     863                        v_ice(ji,jpj-1)  = zFV_prime(ji,jpj-1) / zBV_prime(ji,jpj-1) * vmask(ji,jj,jpj-1) ! last row 
     864                        DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed 
     865                           v_ice(ji,jj)   = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj) / zBV_prime(ji,jj) 
    789866                        END DO             
    790867             
    791868                        !---------------        
    792                         ! 2) Relaxation 
     869                        ! -- Relaxation 
    793870                        !--------------- 
    794                         DO ji = fs_2, fs_jpim1     
    795                            u_ice(ji,jj) = zu_p(ji,jj) + rn_relaxfacU * ( u_ice(ji,jj) - zu_p(ji,jj) ) 
     871                        DO jj = 2, jpj - 1 
     872                           v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) ) 
    796873                        END DO 
    797874 
    798                      END DO ! jj 
     875                     END DO ! ji 
    799876 
    800877                  END DO ! zebra loop 
    801                    
    802                   CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1. ) 
    803                    
    804                ENDIF ! doIterate4u 
    805  
    806                                           ! ---------------------------- ! 
    807                IF ( doIterate4v ) THEN    ! --- Solve for V-velocity --- ! 
    808                                           ! ---------------------------- ! 
    809                                            
    810                   ! what follows could be subroutinized... 
    811                    
    812                   DO k = 1, nn_zebra_step ! "zebra" loop 
    813        
    814                      IF ( k == 1 ); THEN; i_min = 2; ELSE; i_min = 3 
    815           
    816                      zFV(:,:) = 0._wp 
    817  
    818                      DO ji = i_min, fs_jpim1, nn_zebra_step  
    819                          
    820                         !!! HERE WE HAVE ji then jj loops - I don't know what this implies 
    821                         !!! But ZH explain it is critical for convergence speed 
    822  
    823                         !------------------------------------------- 
    824                         ! 1) Tridiagonal system solver for each row 
    825                         !------------------------------------------- 
    826                         ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F 
    827  
    828                         ! --- Right-hand side of tridiagonal system (zFU) 
    829                         DO jj = 2, jpjm1 
    830  
    831                            ! boundary condition substitution (check it is correctly applied !!!) 
    832                            ! see Zhang and Hibler, 1997, Appendix B 
    833                            zAA3 = 0._wp 
    834                            IF ( jj .EQ. 2 )     zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1) 
    835                            IF ( jj .EQ. jpjm1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1) 
    836       
    837                            ! right hand side 
    838                            zFV(ji,jj) = zrhsv(ji,jj) &              ! right-hand side terms 
    839                                &      + zAA3                        ! boundary condition translation 
    840                                &      + zDV(ji,jj) * v_ice(ji-1,jj) ! internal force, j-1 
    841                                &      + zEV(ji,jj) * v_ice(ji+1,jj) ! internal force, j+1 
    842                      
    843                            zFV(ji,jj) = zFV(ji,jj) * vmask(ji,jj)     ! mask - useful ? 
    844          
    845                         END DO 
    846         
    847                         ! --- Thomas Algorithm  
    848                         ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 
    849                         ! Forward sweep 
    850                         DO jj = 3, jpjm1 
    851                            zw         = zAV(ji,jj) / zBV(ji,jj-1) 
    852                            zBV(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 
    853                            zFV(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1)  
    854                         END DO 
    855              
    856                         ! Backward sweep 
    857                         v_ice(ji,jpjm1)  = zFV(ji) / zBV(ji) ! last row 
    858                         DO jj  = jpjm1 - 1, 2, -1 ! can be turned into forward row by substituting jj if needed 
    859                            v_ice(ji,jj)  = zFV(ji) - zCV(ji,jj) * v_ice(ji,jj+1) / zBV(ji,jj) 
    860                         END DO             
    861              
    862                         !---------------        
    863                         ! 2) Relaxation 
    864                         !--------------- 
    865                         DO jj = 2, jpjm1 
    866                            v_ice(ji,jj) = zv_p(ji,jj) + rn_relaxfacV * ( v_ice(ji,jj) - zv_p(ji,jj) ) 
    867                         END DO 
    868  
    869                      END DO ! ji 
    870  
    871                   END DO ! zebra loop 
    872                    
    873                   CALL lbc_lnk( 'icedyn_rhg_vp', v_ice, 'V', -1. )                                           
    874                    
    875                ENDIF ! doIterate4v 
     878                                     
     879               ENDIF !   ll_v_iterate 
     880                
     881               IF ( ( ll_u_iterate .OR. ll_v_iterate ) .OR. jter == nn_nvp ) CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
    876882                               
    877                !-------------------------------------------------- 
    878                ! Continue the loop - based on velocity difference 
    879                !-------------------------------------------------- 
     883               !-------------------------------------------------------------------------------------- 
     884               ! -- Check convergence based on maximum velocity difference, continue or stop the loop 
     885               !-------------------------------------------------------------------------------------- 
    880886 
    881887               !------ 
    882888               ! on U 
    883889               !------ 
    884                IF ( doIterate4u .AND. MOD ( i_inner, nn_lsr_cvgcheck ) == 0 ) THEN 
    885                 
    886                   zuerr(:,:)   = 0._wp    
    887                   DO jj = 2, jpjm1 
    888                      DO ji = 2, jpim1 
    889                         zuerr(ji,jj) = ( u_ice(ji,jj) - zu_p(ji,jj) ) * umask(ji,jj) 
     890               ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check 
     891               ! then we must optimize the use of the mpp_max, which is prohibitive                             
     892                                
     893               IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN 
     894 
     895                  ! - Maximum U-velocity difference                
     896                  zuerr(:,:) = 0._wp 
     897                  DO jj = 2, jpj - 1 
     898                     DO ji = 2, jpi - 1 
     899                        zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) ! * zmask15 ---> MV TEST mask at 15% concentration 
    890900                     END DO 
    891901                  END DO 
    892                    
    893902                  zuerr_max = MAXVAL( zuerr ) 
    894                   CALL mpp_max( 'icedyn_rhg_evp', zuerr_max )   ! max over the global domain 
    895                    
    896                    
    897                   IF ( ln_rhg_chkcvg ) zdiag_u_idif(i_iter) = zuerr_max 
    898                     
    899                   ! "safeguard against bad forcing" -> stop relaxation                     
    900                   IF ( ( i_inner > 1 ) && zuerr_max > rn_uerr_max ) THEN 
    901                       !---> multiply relaxation factor above to stop relaxation 
    902                       zrelaxation_switchU = 0._wp 
     903                  CALL mpp_max( 'icedyn_rhg_evp', zuerr_max )   ! max over the global domain - damned! 
     904                   
     905                  ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 
     906                  IF ( i_inn > 1 && zuerr_max > rn_uerr_max_vp ) THEN 
     907                      IF ( lwp ) " VP rheology error was too large : ", zu_err_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 
     908                      ll_u_iterate = .FALSE. 
    903909                  ENDIF 
    904910                   
    905                   IF ( zuerr_max .LT. rn_uerr_min ) THEN                                         
    906                       doIterate4u = .FALSE. 
    907                       IF ( lwp ) " VP rheology finished in outer U-iteration ", i_outer, " after ", i_inner, " iterations " 
     911                  ! - Stop if error small enough 
     912                  IF ( zuerr_max < rn_uerr_min_vp ) THEN                                         
     913                      IF ( lwp ) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! " 
     914                      ll_u_iterate = .FALSE. 
    908915                  ENDIF 
    909                    
    910                ENDIF 
     916                                                
     917               ENDIF ! ll_u_iterate 
    911918 
    912919               !------ 
     
    914921               !------ 
    915922                
    916                IF ( doIterate4v .AND. MOD ( i_inner, nn_lsr_cvgcheck ) == 0 ) THEN 
    917                 
     923               IF ( ll_v_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN 
     924                
     925                  ! - Maximum V-velocity difference 
    918926                  zverr(:,:)   = 0._wp    
    919                   DO jj = 2, jpjm1 
    920                      DO ji = 2, jpim1 
    921                         zverr(ji,jj) = ( v_ice(ji,jj) - zv_p(ji,jj) ) * vmask(ji,jj) 
     927                  DO jj = 2, jpj - 1 
     928                     DO ji = 2, jpi - 1 
     929                        zverr(ji,jj) = ABS ( ( v_ice(ji,jj) - zv_b(ji,jj) ) ) * vmask(ji,jj,1) 
    922930                     END DO 
    923931                  END DO 
    924932                   
    925933                  zverr_max = MAXVAL( zverr ) 
    926                   CALL mpp_max( 'icedyn_rhg_evp', zverr_max )   ! max over the global domain 
    927                                     
    928                   IF ( ln_rhg_chkcvg ) zdiag_v_idif(i_iter) = zverr_max 
    929                     
    930                   ! "safeguard against bad forcing" -> stop relaxation                     
    931                   IF ( ( i_inner > 1 ) && zverr_max > rn_uerr_max ) THEN 
    932                       !---> multiply relaxation factor above to stop relaxation 
    933                       zrelaxation_switchV = 0._wp 
     934                  CALL mpp_max( 'icedyn_rhg_evp', zverr_max )   ! max over the global domain - damned! 
     935                   
     936                  ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 
     937                  IF ( i_inn > 1 && zverr_max > rn_uerr_max_vp ) THEN 
     938                      IF ( lwp ) " VP rheology error was too large : ", zv_err_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 
     939                      ll_v_iterate = .FALSE. 
    934940                  ENDIF 
    935941                   
    936                   IF ( zverr_max .LT. rn_verr_min ) THEN                                         
    937                       doIterate4v = .FALSE. 
    938                       IF ( lwp ) " VP rheology finished in outer V-iteration ", i_outer, " after ", i_inner, " iterations " 
     942                  ! - Stop if error small enough 
     943                  IF ( zverr_max < rn_verr_min ) THEN                                         
     944                      IF ( lwp ) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! " 
     945                      ll_v_iterate = .FALSE. 
    939946                  ENDIF 
    940947                   
    941                ENDIF 
    942                 
    943                !----------------------------------------------------------------------------- 
    944                ! Convergence diagnostics: residual norm and mean kinetic energy at inner iteration 
    945                !----------------------------------------------------------------------------- 
    946                 
    947                IF( ln_rhg_chkcvg ) THEN 
    948                 
    949                   ! residual norm, velocity difference  
    950                   ! ... zdiag_resnor 
    951                           
    952                   ! mean kinetic energy             
    953                   ztotke = 0._wp 
    954                   DO jj = 2, jpjm1 
    955                      DO ji = 2, jpim1 
    956                         ztotke = ztotke + ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * at_i(ji,jj) & 
    957                            &              ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * z1_32     
    958                      END DO 
    959                   END DO 
    960              
    961                   zdiag_meanke(i_iter) = ztotke / zglob_area 
    962                    
    963                ENDIF 
    964              
    965             ENDIF ! ---    end doIterate4u or doIterate4v 
     948               ENDIF ! ll_v_iterate 
     949                
     950               !--------------------------------------------------------------------------------------- 
     951               ! 
     952               ! --- Calculate extra convergence diagnostics and save them 
     953               ! 
     954               !--------------------------------------------------------------------------------------- 
     955               IF( ln_rhg_chkcvg ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 
     956                          &                         zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 
     957                
     958             
     959            ENDIF ! ---    end ll_u_iterate or ll_v_iterate 
    966960             
    967961         END DO ! i_inn, end of inner loop 
    968962 
    969963      !------------------------------------------------------------------------------! 
     964      ! 
    970965      ! 6) Mask final velocities 
     966      ! 
    971967      !------------------------------------------------------------------------------! 
    972968 
    973969      END ! End of outer loop (i_out) ============================================================================================= 
    974970 
    975       ! 
    976971      !------------------------------------------------------------------------------! 
    977       ! 7) Recompute delta, shear and div (inputs for mechanical redistribution)  
     972      ! 
     973      ! --- Recompute delta, shear and div (inputs for mechanical redistribution)  
     974      ! 
    978975      !------------------------------------------------------------------------------! 
    979       DO jj = 1, jpjm1 
    980          DO ji = 1, jpim1 
     976      ! 
     977      ! OPT: subroutinize ? 
     978             
     979      DO jj = 1, jpj - 1 
     980         DO ji = 1, jpi - 1 
    981981 
    982982            ! shear at F points 
     
    988988      END DO            
    989989       
    990       DO jj = 2, jpjm1 
    991          DO ji = 2, jpim1 ! no vector loop 
     990      DO jj = 2, jpj - 1 
     991         DO ji = 2, jpi - 1 !  
    992992             
    993993            ! tension**2 at T points 
     
    10171017         END DO 
    10181018      END DO 
     1019       
    10191020      CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. ) 
     1021      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    10201022       
    10211023      ! --- Store the stress tensor for the next time step --- ! 
    1022       CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
     1024      ! MV OPT: Are they computed at the end of the tucking fime step ??? 
     1025      ! MV Is stress at previous time step needed for VP, normally no, because they equation is not tucking fime dependent!!! 
     1026      ! 
    10231027      pstress1_i (:,:) = zs1 (:,:) 
    10241028      pstress2_i (:,:) = zs2 (:,:) 
     
    10271031 
    10281032      !------------------------------------------------------------------------------! 
    1029       ! 8) diagnostics 
     1033      ! 
     1034      ! --- Diagnostics 
     1035      ! 
    10301036      !------------------------------------------------------------------------------! 
    1031       ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
     1037      ! MV OPT: subroutinize 
     1038      ! 
     1039      ! --- Ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
    10321040      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 
    10331041         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN 
     
    10461054      ENDIF 
    10471055        
    1048       ! --- divergence, shear and strength --- ! 
     1056      ! --- Divergence, shear and strength --- ! 
    10491057      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
    10501058      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
    10511059      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    10521060 
    1053       ! --- stress tensor --- ! 
     1061      ! --- Stress tensor --- ! 
    10541062      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN 
    10551063         ! 
    10561064         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
    10571065         !          
    1058          DO jj = 2, jpjm1 
    1059             DO ji = 2, jpim1 
     1066         DO jj = 2, jpj - 1 
     1067            DO ji = 2, jpi - 1 
     1068             
    10601069               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    10611070                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  & 
     
    10661075               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    10671076 
    1068 !!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    1069 !!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
    1070 !!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
    1071 !!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    10721077               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    10731078               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
    10741079               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
     1080                
    10751081            END DO 
    10761082         END DO 
     1083 
    10771084         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
    10781085         ! 
     
    10861093 
    10871094         DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
     1095          
    10881096      ENDIF 
    10891097 
    1090       ! --- SIMIP --- ! 
     1098      ! --- SIMIP, terms of tendency for momentum equation --- ! 
    10911099      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
    10921100         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN 
    10931101         ! 
     1102!!!!!!!!!         ATTENTION LA FORCE INTERNE DOIT ETTRE RECALCIULEEE ICCI !!!!!!!!!!!!!!!! 
    10941103         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., & 
    10951104            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. ) 
     
    11091118            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) ) 
    11101119         ! 
    1111          DO jj = 2, jpjm1 
    1112             DO ji = 2, jpim1 
     1120         DO jj = 2, jpj - 1 
     1121            DO ji = 2, jpi - 1 
    11131122               ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    11141123               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj) 
     
    11431152      ENDIF 
    11441153 
    1145       ! 
    1146       ! --- convergence diagnostics --- ! 
     1154      ! --- Convergence diagnostics --- ! 
    11471155      IF( ln_rhg_chkcvg ) THEN 
    11481156           
    1149           ! still miss residual norm 
    1150                 
    1151          IF ( iom_use('mean_tke_ice') ) THEN 
    1152             CALL iom_put( 'mean_tke_ice', zdiag_meanke )                   
    1153             DEALLOCATE( zdiag_meanke )        
    1154          ENDIF 
    1155           
    1156          IF ( iom_use('u_veldif') ) THEN 
    1157             CALL iom_put ( zdiag_u_idif(i_iter) ) 
    1158             DEALLOCATE( zdiag_u_idif ) 
    1159          ENDIF 
    1160           
    1161          IF ( iom_use('v_veldif') ) THEN 
    1162             CALL iom_put ( zdiag_v_idif(i_iter) ) 
    1163             DEALLOCATE( zdiag_v_idif ) 
    1164          ENDIF 
    1165           
     1157         IF( iom_use('uice_cvg') ) THEN 
     1158            CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , & 
     1159                  &                        ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     1160         ENDIF    
     1161         
    11661162      ENDIF ! ln_rhg_chkcvg 
    11671163 
    1168 ! EVP convergence diag       
    1169 !         IF( iom_use('uice_cvg') ) THEN 
    1170 !          
    1171 !         
    1172 !           ! EVP-diag 
    1173 !           IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
    1174 !              CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 
    1175 !                 &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 
    1176 !           ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
    1177 !              CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 
    1178 !                 &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
    1179 !           ENDIF 
    1180  
    1181          ENDIF 
    1182       ENDIF       
    11831164      ! 
    11841165      DEALLOCATE( zmsk00, zmsk15 ) 
    1185  
    11861166 
    11871167   END SUBROUTINE ice_dyn_rhg_vp 
    11881168    
     1169    
     1170    
     1171   SUBROUTINE rhg_cvg_vp( kt, kiter, kitermax, pu, pv, pmt, puerr_max, pverr_max, pglob_area, & 
     1172                  &       prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV ) 
     1173    
     1174!                  CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 
     1175!                          &        zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 
     1176      !!---------------------------------------------------------------------- 
     1177      !!                    ***  ROUTINE rhg_cvg_vp  *** 
     1178      !!                      
     1179      !! ** Purpose :   check convergence of VP ice rheology 
     1180      !! 
     1181      !! ** Method  :   create a file ice_cvg.nc containing a few convergence diagnostics 
     1182      !!                This routine is called every sub-iteration, so it is cpu expensive 
     1183      !! 
     1184      !!                Calculates / stores 
     1185      !!                   - maximum absolute U-V difference (uice_cvg, u_dif, v_dif, m/s) 
     1186      !!                   - residuals in U, V and UV-mean taken as square-root of area-weighted mean square residual (u_res, v_res, vel_res, N/m2) 
     1187      !!                   - mean kinetic energy (mke_ice, J/m2) 
     1188      !!                    
     1189      !! 
     1190      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)    
     1191      !!---------------------------------------------------------------------- 
     1192      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax      ! ocean time-step index 
     1193      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pmt              ! now velocity and mass per unit area  
     1194      REAL(wp),                 INTENT(in) ::   puerr_max, pverr_max     ! absolute mean velocity difference 
     1195      REAL(wp),                 INTENT(in) ::   pglob_area               ! global ice area 
     1196      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsu, pAU, pBU, pCU, pDU, pEU ! linear system coefficients  
     1197      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsv, pAV, pBV, pCV, pDV, pEV 
     1198      !! 
     1199      INTEGER           ::   it, idtime, istatus 
     1200      INTEGER           ::   ji, jj          ! dummy loop indices 
     1201      REAL(wp)          ::   zveldif, zu_res_mean, zv_res_mean, zmke, zu, zv ! local scalars 
     1202      REAL(wp), DIMENSION(jpi,jpj)   ::   zu_res(:,:), zv_res(:,:), zvel2(:,:) ! local arrays 
     1203                                                                              
     1204      CHARACTER(len=20) ::   clname 
     1205      ::   zres           ! check convergence 
     1206      !!---------------------------------------------------------------------- 
     1207 
     1208      ! create file 
     1209      IF( kt == nit000 .AND. kiter == 1 ) THEN 
     1210         ! 
     1211         IF( lwp ) THEN 
     1212            WRITE(numout,*) 
     1213            WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 
     1214            WRITE(numout,*) '~~~~~~~~~~~' 
     1215         ENDIF 
     1216         ! 
     1217         IF( lwm ) THEN 
     1218            clname = 'ice_cvg.nc' 
     1219            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 
     1220            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 
     1221            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime ) 
     1222            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid_ucvg ) 
     1223            istatus = NF90_DEF_VAR( ncvgid, 'u_res', NF90_DOUBLE   , (/ idtime /), nvarid_ures ) 
     1224            istatus = NF90_DEF_VAR( ncvgid, 'v_res', NF90_DOUBLE   , (/ idtime /), nvarid_vres ) 
     1225            istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE   , (/ idtime /), nvarid_velres ) 
     1226            istatus = NF90_DEF_VAR( ncvgid, 'u_dif', NF90_DOUBLE   , (/ idtime /), nvarid_udif ) 
     1227            istatus = NF90_DEF_VAR( ncvgid, 'v_dif', NF90_DOUBLE   , (/ idtime /), nvarid_vdif ) 
     1228            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE   , (/ idtime /), nvarid_mke ) 
     1229            istatus = NF90_ENDDEF(ncvgid) 
     1230         ENDIF 
     1231         ! 
     1232      ENDIF 
     1233 
     1234      ! time 
     1235      it = ( kt - 1 ) * kitermax + kiter 
     1236 
     1237      ! --- Max absolute velocity difference with previous iterate (zveldif) 
     1238      ! EVP code       
     1239!      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 
     1240!         zveldif = 0._wp 
     1241!      ELSE 
     1242!         DO jj = 1, jpj 
     1243!            DO ji = 1, jpi 
     1244!               zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     1245!                  &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 
     1246!            END DO 
     1247!         END DO 
     1248!         zveldif = MAXVAL( zres ) 
     1249!         CALL mpp_max( 'icedyn_rhg_vp', zveldif )   ! max over the global domain 
     1250!      ENDIF 
     1251      ! VP code 
     1252      zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code  
     1253                                            ! if puerrmask and pverrmax are masked at 15% (TEST) 
     1254       
     1255      ! -- Mean residual (N/m^2), zu_res_mean 
     1256      ! Here we take the residual of the linear system (N/m^2),  
     1257      ! We define it as in mitgcm: square-root of area-weighted mean square residual 
     1258      ! Local residual r = Ax - B expresses to which extent the momentum balance is verified  
     1259      ! i.e., how close we are to a solution 
     1260      DO jj = 2, jpj - 1 
     1261         DO ji = 2, jpi - 1                                       
     1262            zu_res(ji,jj)  = zrhsu(ji,jj) + zDU(ji,jj) * pu(ji,jj-1) + zEU(ji,jj) * pu(ji,jj+1) & 
     1263               &           - zAU(ji,jj) * pu(ji-1,jj) - zBU(ji,jj) * pu(ji,jj) - zCU(ji,jj) * pu(ji+1,jj) 
     1264                            
     1265            zv_res(ji,jj)  = zrhsv(ji,jj) + zDV(ji,jj) * pv(ji-1,jj) + zEV(ji,jj) * pv(ji+1,jj) & 
     1266               &           - zAV(ji,jj) * pv(ji,jj-1) - zBV(ji,jj) * pv(ji,jj) - zCV(ji,jj) * pv(ji,jj+1)                            
     1267          END DO 
     1268       END DO                   
     1269       zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:) ) 
     1270       zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:) ) 
     1271       zu_res_mean = SQRT( zu_resmean / pglob_area ) 
     1272       zv_res_mean = SQRT( zv_resmean / pglob_area ) 
     1273       zvelres     = MEAN( zu_res_mean, zv_res_mean ) 
     1274                                           
     1275       ! -- Global mean kinetic energy per unit area (J/m2)   
     1276       DO jj = 2, jpj - 1 
     1277          DO ji = 2, jpi - 1                    
     1278             zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 
     1279             zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 
     1280             zvel2(:,:)  = zu*zu + zv*zv              ! square of ice velocity at T-point   
     1281          END DO 
     1282       END DO 
     1283        
     1284       zmke = 0.5_wp * globsum( 'ice_rhg_vp', zmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) 
     1285                   
     1286         !                                                ! ==================== ! 
     1287 
     1288      IF( lwm ) THEN 
     1289         ! write variables 
     1290            istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) ) 
     1291            istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) 
     1292            istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) 
     1293            istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) 
     1294            istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) 
     1295            istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) 
     1296            istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) 
     1297         ! close file 
     1298         IF( kt == nitend )   istatus = NF90_CLOSE( ncvgid ) 
     1299      ENDIF 
     1300       
     1301   END SUBROUTINE rhg_cvg_vp 
     1302    
     1303 
    11891304 
    11901305   SUBROUTINE rhg_vp_rst( cdrw, kt ) 
Note: See TracChangeset for help on using the changeset viewer.