New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_rhg_vp.F90 in NEMO/branches/2020/SI3-03_VP_rheology/src/ICE – NEMO

source: NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90 @ 15527

Last change on this file since 15527 was 15527, checked in by vancop, 7 months ago

Optimize comms for VP

File size: 86.9 KB
Line 
1MODULE icedyn_rhg_vp
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_vp  ***
4   !!   Sea-Ice dynamics : Viscous-plastic rheology with LSR technique
5   !!======================================================================
6   !!
7   !! History :   -   !  1997-20  (J. Zhang, M. Losch) Original code, implementation into mitGCM
8   !!            4.0  !  2020-09  (M. Vancoppenolle) Adaptation to SI3
9   !!
10   !! CLEAN-UP 2 — remove iomputs
11   !!----------------------------------------------------------------------
12#if defined key_si3
13   !!----------------------------------------------------------------------
14   !!   'key_si3'                                       SI3 sea-ice model
15   !!----------------------------------------------------------------------
16   !!   ice_dyn_rhg_vp : computes ice velocities from VP rheolog with LSR solvery
17   !!----------------------------------------------------------------------
18   USE phycst         ! Physical constants
19   USE dom_oce        ! Ocean domain
20   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
21   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
22   USE ice            ! sea-ice: ice variables
23   USE icevar         ! ice_var_sshdyn
24   USE icedyn_rdgrft  ! sea-ice: ice strength
25   USE bdy_oce , ONLY : ln_bdy 
26   USE bdyice 
27#if defined key_agrif
28   USE agrif_ice_interp
29#endif
30   !
31   USE in_out_manager ! I/O manager
32   USE iom            ! I/O manager library
33   USE lib_mpp        ! MPP library
34   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
35   USE lbclnk         ! lateral boundary conditions (or mpp links)
36   USE prtctl         ! Print control
37
38   USE netcdf         ! NetCDF library for convergence test
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   ice_dyn_rhg_vp   ! called by icedyn_rhg.F90
43
44   !!
45   INTEGER ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp)     
46
47   !! for convergence tests
48   INTEGER ::   ncvgid        ! netcdf file id
49   INTEGER ::   nvarid_ures, nvarid_vres, nvarid_velres
50   INTEGER ::   nvarid_uerr_max, nvarid_verr_max, nvarid_velerr_max
51   INTEGER ::   nvarid_umad, nvarid_vmad, nvarid_velmad
52   INTEGER ::   nvarid_umad_outer, nvarid_vmad_outer, nvarid_velmad_outer
53   INTEGER ::   nvarid_mke
54
55   !INTEGER ::   nvarid_ures_xy, nvarid_vres_xy
56
57         !!! -> store velocity at last inner iteration zu_end_outer, zv_end_outer
58         !!! -> send it to rhg_cvg
59         !!! -> for 2nd and larger iteration, calculate mad_outer and store it (rhg_cvg)
60         !!! -> store spatially dependent residual (rhg_cvg)
61         !!! -> store spatially dependent mad_outer (rhg_cvg)
62
63   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
64
65   !!----------------------------------------------------------------------
66   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
67   !! $Id: icedyn_rhg_vp.F90 13279 2020-07-09 10:39:43Z clem $
68   !! Software governed by the CeCILL license (see ./LICENSE)
69   !!----------------------------------------------------------------------
70   
71CONTAINS
72
73   SUBROUTINE ice_dyn_rhg_vp( kt, pshear_i, pdivu_i, pdelta_i )
74      !!-------------------------------------------------------------------
75      !!
76      !!                 ***  SUBROUTINE ice_dyn_rhg_vp  ***
77      !!                             VP-LSR-C-grid
78      !!
79      !! ** Purpose : determines sea ice drift from wind stress, ice-ocean
80      !!  stress and sea-surface slope. Internal forces assume viscous-plastic rheology (Hibler, 1979)
81      !!
82      !! ** Method
83      !! 
84      !!  The resolution algorithm  follows from Zhang and Hibler (1997) and Losch (2010)
85      !!  with elements from Lemieux and Tremblay (2008) and Lemieux and Tremblay (2009)
86      !! 
87      !!  The components of the momentum equations are arranged following the ideas of Zhang and Hibler (1997)
88      !!
89      !!  f1(u) = g1(v)
90      !!  f2(v) = g2(u)
91      !!
92      !!  The right-hand side (RHS) is explicit
93      !!  The left-hand side (LHS) is implicit
94      !!  Coriolis is part of explicit terms, whereas ice-ocean drag is implicit
95      !!
96      !!  Two iteration levels (outer and inner loops) are used to solve the equations
97      !!
98      !!  The outer loop (OL, typically 10 iterations) is there to deal with the (strong) non-linearities in the equation
99      !!
100      !!  The inner loop (IL, typically 1500 iterations) is there to solve the linear problem with a line-successive-relaxation algorithm
101      !!
102      !!  The velocity used in the non-linear terms uses a "modified euler time step" (not sure its the correct term),
103      !!! with uk = ( uk-1 + uk-2 ) / 2.
104      !! 
105      !!  * Spatial discretization
106      !!
107      !!  Assumes a C-grid
108      !!
109      !!  The points in the C-grid look like this, my darling
110      !!
111      !!                              (ji,jj)
112      !!                                 |
113      !!                                 |
114      !!                      (ji-1,jj)  |  (ji,jj)
115      !!                             ---------   
116      !!                            |         |
117      !!                            | (ji,jj) |------(ji,jj)
118      !!                            |         |
119      !!                             ---------   
120      !!                     (ji-1,jj-1)     (ji,jj-1)
121      !!
122      !! ** Inputs  : - wind forcing (stress), oceanic currents
123      !!                ice total volume (vt_i) per unit area
124      !!                snow total volume (vt_s) per unit area
125      !!
126      !! ** Action  :
127      !!             
128      !! ** Steps   :
129      !!           
130      !! ** Notes   :
131      !!             
132      !! References : Zhang and Hibler, JGR 1997; Losch et al., OM 2010., Lemieux et al., 2008, 2009, ... 
133      !!             
134      !!             
135      !!-------------------------------------------------------------------
136      !!
137      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
138      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
139      !!
140      LOGICAL ::   ll_u_iterate, ll_v_iterate   ! continue iteration or not
141      !
142      INTEGER ::   ji, ji2, jj, jj2, jn          ! dummy loop indices
143      INTEGER ::   i_out, i_inn, i_inn_tot  !
144      INTEGER ::   ji_min, jj_min      !
145      INTEGER ::   nn_zebra_vp         ! number of zebra steps
146
147      !
148      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
149      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
150      REAL(wp) ::   zglob_area                                          ! global ice area for diagnostics
151      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
152      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass and volume
153      REAL(wp) ::   zds2, zdt, zdt2, zdiv, zdiv2                        ! temporary scalars
154      REAL(wp) ::   zp_delstar_f                                      !
155      REAL(wp) ::   zu_cV, zv_cU                                        !
156      REAL(wp) ::   zfac, zfac1, zfac2, zfac3
157      REAL(wp) ::   zt12U, zt11U, zt22U, zt21U, zt122U, zt121U
158      REAL(wp) ::   zt12V, zt11V, zt22V, zt21V, zt122V, zt121V
159      REAL(wp) ::   zAA3, zw, ztau, zuerr_max, zverr_max
160      !
161      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice
162      REAL(wp), DIMENSION(jpi,jpj) ::   za_iU  , za_iV                      ! ice fraction on U/V points
163      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! Acceleration term contribution to RHS
164      REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step
165      !
166      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltat, zdelstar_t           ! Delta & Delta* at T-points
167      REAL(wp), DIMENSION(jpi,jpj) ::   ztens, zshear                   ! Tension, shear
168      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delstar_t                  ! P/delta* at T points
169      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zet                        ! Viscosity pre-factors at T points
170      REAL(wp), DIMENSION(jpi,jpj) ::   zef                             ! Viscosity pre-factor at F point
171      !
172      REAL(wp), DIMENSION(jpi,jpj) ::   zmt                             ! Mass per unit area at t-point
173      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! Coriolis factor (m*f) at t-point
174      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points
175      REAL(wp), DIMENSION(jpi,jpj) ::   zu_c, zv_c                      ! "current" ice velocity (m/s), average of previous two OL iterates
176      REAL(wp), DIMENSION(jpi,jpj) ::   zu_b, zv_b                      ! velocity at previous sub-iterate
177      REAL(wp), DIMENSION(jpi,jpj) ::   zuerr, zverr                    ! absolute U/Vvelocity difference between current and previous sub-iterates
178      !
179      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
180      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
181      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
182      !                                                                 !    ice bottom surface if ice is embedded   
183      REAL(wp), DIMENSION(jpi,jpj) ::   zCwU, zCwV                      ! ice-ocean drag pre-factor (rho*c*module(u))
184      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
185      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
186      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
187      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi_rhsu, ztauy_oi_rhsv    ! ice-ocean stress RHS contribution at U-V points
188      REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsu, zs2_rhsu, zs12_rhsu   ! internal stress contributions to RHSU
189      REAL(wp), DIMENSION(jpi,jpj) ::   zs1_rhsv, zs2_rhsv, zs12_rhsv   ! internal stress contributions to RHSV
190      REAL(wp), DIMENSION(jpi,jpj) ::   zf_rhsu, zf_rhsv                ! U- and V- components of internal force RHS contributions
191      REAL(wp), DIMENSION(jpi,jpj) ::   zrhsu, zrhsv                    ! U and V RHS
192      REAL(wp), DIMENSION(jpi,jpj) ::   zAU, zBU, zCU, zDU, zEU         ! Linear system coefficients, U equation
193      REAL(wp), DIMENSION(jpi,jpj) ::   zAV, zBV, zCV, zDV, zEV         ! Linear system coefficients, V equation
194      REAL(wp), DIMENSION(jpi,jpj) ::   zFU, zFU_prime, zBU_prime       ! Rearranged linear system coefficients, U equation
195      REAL(wp), DIMENSION(jpi,jpj) ::   zFV, zFV_prime, zBV_prime       ! Rearranged linear system coefficients, V equation
196!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
197!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
198     !
199      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! mask for lots of ice (1), little ice (0)
200      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence (1), no ice (0)
201      !
202      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
203      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small
204      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
205      !! --- diags
206      REAL(wp)                     ::   zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y
207      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f           ! stress tensor components
208      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
209      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztaux_oi, ztauy_oi
210      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice, zdiag_ymtrp_ice ! X/Y-component of ice mass transport (kg/s, SIMIP)
211      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw, zdiag_ymtrp_snw ! X/Y-component of snow mass transport (kg/s, SIMIP)
212      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp, zdiag_yatrp         ! X/Y-component of area transport (m2/s, SIMIP)
213
214      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zvel_res                         ! Residual of the linear system at last iteration
215      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zvel_diff                        ! Absolute velocity difference @last outer iteration
216                       
217      !!----------------------------------------------------------------------------------------------------------------------
218      ! sea surface height
219      !    embedded sea ice: compute representative ice top surface
220      !    non-embedded sea ice: use ocean surface for slope calculation
221      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
222
223      !!----------------------------------------------------------------------------------------------------------------------
224
225      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)'
226      IF( lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)'
227           
228      !------------------------------------------------------------------------------!
229      !
230      ! --- Initialization
231      !
232      !------------------------------------------------------------------------------!
233     
234      ! for diagnostics and convergence tests
235      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) )
236      DO jj = 1, jpj
237         DO ji = 1, jpi
238            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
239            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
240         END DO
241      END DO
242     
243      IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2
244                         ELSE; nn_zebra_vp = 1; ENDIF
245
246      nn_nvp = nn_nout_vp * nn_ninn_vp ! maximum number of iterations
247
248      IF( lwp )   WRITE(numout,*) ' ln_zebra_vp : ', ln_zebra_vp
249      IF( lwp )   WRITE(numout,*) ' nn_zebra_vp : ', nn_zebra_vp
250      IF( lwp )   WRITE(numout,*) ' nn_nvp      : ', nn_nvp
251     
252      zrhoco = rau0 * rn_cio 
253
254      ! ecc2: square of yield ellipse eccentricity
255      ecc2    = rn_ecc * rn_ecc
256      z1_ecc2 = 1._wp / ecc2
257               
258      ! Initialise convergence checks
259      IF( ln_rhg_chkcvg ) THEN
260
261         ! allocate spatially-resolved convergence diagnostics
262         ALLOCATE( zvel_res(jpi,jpj), zvel_diff(jpi,jpj) )
263     
264         ! ice area for proper weighted mean averages
265         zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! global ice area (km2)
266         
267      ENDIF
268
269      ! Landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
270      ! MV: Not working yet...
271      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile
272      ELSE                         ;   zkt = 0._wp
273      ENDIF
274
275      zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp
276      zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp
277      zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp;
278      zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp;
279      zrhsu  (:,:) = 0._wp ; zrhsv  (:,:) = 0._wp
280      zf_rhsu(:,:) = 0._wp ; zf_rhsv(:,:) = 0._wp
281
282      !------------------------------------------------------------------------------!
283      !
284      ! --- Time-independent quantities
285      !
286      !------------------------------------------------------------------------------!
287     
288      CALL ice_strength ! strength at T points
289     
290      !---------------------------
291      ! -- F-mask (code from EVP)
292      !---------------------------
293      ! MartinV:
294      ! In EVP routine, zfmask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip)
295      ! I am not sure the same recipe applies here
296     
297      ! - ocean/land mask
298      DO jj = 1, jpj - 1
299         DO ji = 1, jpi - 1
300            zfmask(ji,jj)    = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
301         END DO
302      END DO
303
304      ! Lateral boundary conditions on velocity (modify zfmask)
305      ! Can be computed once for all, at first time step, for all rheologies
306      DO jj = 2, jpj - 1
307         DO ji = 2, jpi - 1   ! vector opt.
308            IF( zfmask(ji,jj) == 0._wp ) THEN
309               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
310                  &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
311            ENDIF
312         END DO
313      END DO
314      DO jj = 2, jpj - 1
315         IF( zfmask(1,jj) == 0._wp ) THEN
316            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) )
317         ENDIF
318         IF( zfmask(jpi,jj) == 0._wp ) THEN
319            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpi - 1,jj,1), umask(jpi,jj-1,1) ) )
320        ENDIF
321      END DO
322      DO ji = 2, jpi - 1
323         IF( zfmask(ji,1) == 0._wp ) THEN
324            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) )
325         ENDIF
326         IF( zfmask(ji,jpj) == 0._wp ) THEN
327            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpj - 1,1) ) )
328         ENDIF
329      END DO
330     
331      CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp )
332     
333      !----------------------------------------------------------------------------------------------------------
334      ! -- Time-independent pre-factors for acceleration, ocean drag, coriolis, atmospheric drag, surface tilt
335      !----------------------------------------------------------------------------------------------------------
336      ! Compute all terms & factors independent of velocities, or only depending on velocities at previous time step     
337
338      zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:)       ! Snow and ice mass at T-point
339      zmf(:,:) = zmt(:,:) * ff_t(:,:)                      ! Coriolis factor at T points (m*f)
340     
341      DO jj = 2, jpj - 1
342         DO ji = 2, jpi - 1
343
344            ! Ice fraction at U-V points
345            za_iU(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
346            za_iV(ji,jj)    = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
347
348            ! Snow and ice mass at U-V points
349            zm1             = zmt(ji,jj)
350            zm2             = zmt(ji+1,jj)
351            zm3             = zmt(ji,jj+1)
352            zmassU          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
353            zmassV          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
354                         
355            ! Mass per unit area divided by time step
356            zmassU_t(ji,jj) = zmassU * r1_rdtice
357            zmassV_t(ji,jj) = zmassV * r1_rdtice
358
359            ! Acceleration term contribution to RHS (depends on velocity at previous time step)           
360            zmU_t(ji,jj)    = zmassU_t(ji,jj) * u_ice(ji,jj)
361            zmV_t(ji,jj)    = zmassV_t(ji,jj) * v_ice(ji,jj)
362           
363            ! Ocean currents at U-V points
364            v_oceU(ji,jj)   = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
365            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
366           
367            ! Wind stress
368            ztaux_ai(ji,jj) = za_iU(ji,jj) * utau_ice(ji,jj)
369            ztauy_ai(ji,jj) = za_iV(ji,jj) * vtau_ice(ji,jj)
370
371            ! Force due to sea surface tilt(- m*g*GRAD(ssh))
372            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
373            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
374
375            ! Mask for ice presence (1) or absence (0)
376            zmsk00x(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
377            zmsk00y(ji,jj)  = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
378
379            ! Mask for lots of ice (1) or little ice (0)
380            IF ( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
381            ELSE                                                      ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
382            IF ( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
383            ELSE                                                      ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF             
384
385         END DO
386      END DO   
387           
388      !------------------------------------------------------------------------------!
389      !
390      ! --- Start outer loop
391      !
392      !------------------------------------------------------------------------------!
393     
394      zu_c(:,:) = u_ice(:,:)
395      zv_c(:,:) = v_ice(:,:)
396     
397      i_inn_tot = 0
398
399      DO i_out = 1, nn_nout_vp
400
401
402         ! Velocities used in the non linear terms are the average of the past two iterates
403         ! u_it = 0.5 * ( u_{it-1} + u_{it-2} )
404         ! Also used in Hibler and Ackley (1983); Zhang and Hibler (1997); Lemieux and Tremblay (2009)
405         zu_c(:,:)       = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) )
406         zv_c(:,:)       = 0.5_wp * ( v_ice(:,:) + zv_c(:,:) )
407         
408         !------------------------------------------------------------------------------!
409         !
410         ! --- Right-hand side (RHS) of the linear problem
411         !
412         !------------------------------------------------------------------------------!
413         ! In the outer loop, one needs to update all RHS terms
414         ! with explicit velocity dependencies (viscosities, coriolis, ocean stress)
415         ! as a function of "current" velocities (uc, vc)
416     
417         !------------------------------------------
418         ! -- Strain rates, viscosities and P/Delta
419         !------------------------------------------
420
421         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
422         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
423            DO ji = 1, jpi - 1
424
425               ! shear at F points
426               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)   &
427                  &         + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
428                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
429
430            END DO
431         END DO
432
433         CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! necessary, zds2 uses jpi/jpj values for zds
434
435         DO jj = 2, jpj        ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
436            DO ji = 2, jpi     !
437
438               ! shear**2 at T points (doc eq. A16)
439               zds2  = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
440                  &    + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
441                  &    ) * 0.25_wp * r1_e1e2t(ji,jj)
442             
443               ! divergence at T points
444               zdiv  = ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj)   &
445                  &    + e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1)   &
446                  &    ) * r1_e1e2t(ji,jj)
447               zdiv2 = zdiv * zdiv
448               
449               ! tension at T points
450               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)   &
451                  &    - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
452                  &    ) * r1_e1e2t(ji,jj)
453               zdt2  = zdt * zdt
454               
455               ! delta at T points
456               zdeltat(ji,jj)        = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
457               
458               ! delta* at T points (following Lemieux and Dupont, GMD 2020)
459               zdelstar_t(ji,jj)     = zdeltat(ji,jj) + rn_creepl ! OPT zdelstar_t can be totally removed and put into next line directly. Could change results
460             
461               ! P/delta* at T-points
462               zp_delstar_t(ji,jj)   = strength(ji,jj) / zdelstar_t(ji,jj)
463               
464               ! Temporary zzt and zet factors at T-points
465               zzt(ji,jj)            = zp_delstar_t(ji,jj) * r1_e1e2t(ji,jj)
466               zet(ji,jj)            = zzt(ji,jj)     * z1_ecc2 
467                         
468            END DO
469         END DO
470         
471         CALL lbc_lnk( 'icedyn_rhg_vp', zp_delstar_t , 'T', 1. ) ! necessary, used for ji = 1 and jj = 1
472
473         DO jj = 1, jpj - 1
474            DO ji = 1, jpi - 1
475               
476               ! P/delta* at F points
477               zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) )
478               
479               ! Temporary zef factor at F-point
480               zef(ji,jj)      = zp_delstar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) * 0.5_wp 
481
482            END DO
483         END DO
484
485         !---------------------------------------------------
486         ! -- Ocean-ice drag and Coriolis RHS contributions
487         !---------------------------------------------------
488
489         DO jj = 2, jpj - 1
490             DO ji = 2, jpi - 1
491               
492                !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
493                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)
494                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)
495               
496                !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
497                zCwU(ji,jj)          = za_iU(ji,jj) * zrhoco * SQRT( ( zu_c (ji,jj) - u_oce (ji,jj) ) * ( zu_c (ji,jj) - u_oce (ji,jj) )  &
498                  &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
499                zCwV(ji,jj)          = za_iV(ji,jj) * zrhoco * SQRT( ( zv_c (ji,jj) - v_oce (ji,jj) ) * ( zv_c (ji,jj) - v_oce (ji,jj) )  &
500                  &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
501                 
502                !--- Ocean-ice drag contributions to RHS
503                ztaux_oi_rhsu(ji,jj) = zCwU(ji,jj) * u_oce(ji,jj)
504                ztauy_oi_rhsv(ji,jj) = zCwV(ji,jj) * v_oce(ji,jj)
505               
506                ! --- U-component of Coriolis Force (energy conserving formulation)
507                zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  &
508                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * zv_c(ji  ,jj) + e1v(ji  ,jj-1) * zv_c(ji  ,jj-1) )  &
509                           &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_c(ji+1,jj) + e1v(ji+1,jj-1) * zv_c(ji+1,jj-1) ) )
510                           
511                ! --- V-component of Coriolis Force (energy conserving formulation)
512                zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  &
513                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * zu_c(ji,jj  ) + e2u(ji-1,jj  ) * zu_c(ji-1,jj  ) )  &
514                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) )
515
516             END DO
517         END DO
518
519         !-------------------------------------
520         ! -- Internal stress RHS contribution
521         !-------------------------------------
522
523         ! --- Stress contributions at T-points         
524         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1 & zs2
525            DO ji = 2, jpi !
526           
527               ! sig1 contribution to RHS of U-equation at T-points
528               zs1_rhsu(ji,jj) =   zzt(ji,jj) * ( e1v(ji,jj)    * zv_c(ji,jj) - e1v(ji,jj-1)    * zv_c(ji,jj-1) )   &
529                               &                - zp_delstar_t(ji,jj) * zdeltat(ji,jj)
530                                           
531               ! sig2 contribution to RHS of U-equation at T-points           
532               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) 
533
534               ! sig1 contribution to RHS of V-equation at T-points
535               zs1_rhsv(ji,jj) =   zzt(ji,jj) * ( e2u(ji,jj)    * zu_c(ji,jj) - e2u(ji-1,jj)    * zu_c(ji-1,jj) )   & 
536                               &                - zp_delstar_t(ji,jj) * zdeltat(ji,jj)
537
538               ! sig2 contribution to RHS of V-equation  at T-points
539               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)
540
541            END DO
542         END DO
543         
544         ! --- Stress contributions at F-points         
545         ! MV NOTE: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing
546         ! My guess is that this is the way to enforce boundary conditions on strain rate tensor
547
548         DO jj = 1, jpj - 1
549            DO ji = 1, jpi - 1
550               
551               ! sig12 contribution to RHS of U equation at F-points
552               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)
553               
554               ! sig12 contribution to RHS of V equation at F-points
555               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)
556
557            END DO
558         END DO
559
560         ! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002)
561         ! OPT: merge with next loop and use intermediate scalars for zf_rhsu
562         
563         DO jj = 2, jpj - 1
564            DO ji = 2, jpi - 1               
565
566               ! --- U component of internal force contribution to RHS at U points
567               zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * & 
568                              (    e2u(ji,jj)    * ( zs1_rhsu(ji+1,jj) - zs1_rhsu(ji,jj) )                                                                 &
569                  &      +         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) )     &
570                  &      + 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) ) )
571                 
572               ! --- V component of internal force contribution to RHS at V points
573               zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * &
574                  &           (    e1v(ji,jj)    * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) )                                                                 &
575                  &      -         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) )     &
576                  &      + 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) ) )
577
578            END DO
579         END DO
580         
581         !---------------------------
582         ! -- Sum RHS contributions
583         !---------------------------
584         !
585         ! OPT: could use intermediate scalars to reduce memory access
586         DO jj = 2, jpj - 1
587            DO ji = 2, jpi - 1
588           
589               zrhsu(ji,jj) = zmU_t(ji,jj) + ztaux_ai(ji,jj) + ztaux_oi_rhsu(ji,jj) + zspgU(ji,jj) + zCorU(ji,jj) + zf_rhsu(ji,jj)
590               zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsv(ji,jj)
591
592            END DO
593         END DO
594         
595         !---------------------------------------------------------------------------------------!
596         !
597         ! --- Linear system matrix
598         !
599         !---------------------------------------------------------------------------------------!
600     
601         ! Linear system matrix contains all implicit contributions
602         ! 1) internal forces, 2) acceleration, 3) ice-ocean drag
603
604         ! The linear system equation is written as follows
605         ! AU * u_{i-1,j} + BU * u_{i,j}   + CU * u_{i+1,j}
606         !                = DU * u_{i,j-1} + EU * u_{i,j+1} + RHS          (! my convention, not the same as ZH97 )         
607         
608         ! MV Note 1: martin losch applies boundary condition to BU in mitGCM - check whether it is necessary here ?
609         ! MV Note 2: "T" factor calculations can be optimized by putting things out of the loop
610         !         only zzt and zet are iteration-dependent, other only depend on scale factors
611                 
612         DO ji = 2, jpi - 1 ! internal domain do loop
613            DO jj = 2, jpj - 1
614
615               !-------------------------------------
616               ! -- Internal forces LHS contribution
617               !-------------------------------------
618               !
619               ! --- U-component
620               !
621               ! "T" factors (intermediate results)
622               !
623               zfac       = 0.5_wp * r1_e1e2u(ji,jj)
624               zfac1      =         zfac * e2u(ji,jj)
625               zfac2      =         zfac * r1_e2u(ji,jj)
626               zfac3      = 2._wp * zfac * r1_e1u(ji,jj)
627
628               zt11U      =   zfac1 * zzt(ji,jj)
629               zt12U      =   zfac1 * zzt(ji+1,jj)
630         
631               zt21U      =   zfac2 * zet(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)
632               zt22U      =   zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj)
633         
634               zt121U     =   zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)
635               zt122U     =   zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)
636               
637               !
638               ! Linear system coefficients
639               !
640               zAU(ji,jj) = -   zt11U           * e2u(ji-1,jj) -   zt21U           * r1_e2u(ji-1,jj)
641               zBU(ji,jj) =   ( zt11U + zt12U ) * e2u(ji,jj)   + ( zt21U + zt22U ) * r1_e2u(ji,jj)   + ( zt121U + zt122U ) * r1_e1u(ji,jj)
642               zCU(ji,jj) = -   zt12U           * e2u(ji+1,jj) -   zt22U           * r1_e2u(ji+1,jj)
643         
644               zDU(ji,jj) =     zt121U * r1_e1u(ji,jj-1)
645               zEU(ji,jj) =     zt122U * r1_e1u(ji,jj+1)
646             
647               !
648               ! --- V-component
649               !
650               ! "T" factors (intermediate results)
651               !
652               zfac       = 0.5_wp * r1_e1e2v(ji,jj)
653               zfac1      =         zfac * e1v(ji,jj)
654               zfac2      =         zfac * r1_e1v(ji,jj)
655               zfac3      = 2._wp * zfac * r1_e2v(ji,jj)
656         
657               zt11V      =   zfac1 * zzt(ji,jj)
658               zt12V      =   zfac1 * zzt(ji,jj+1)
659         
660               zt21V      =   zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)
661               zt22V      =   zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1)
662         
663               zt121V     =   zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)
664               zt122V     =   zfac3 * zef(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)
665
666               !
667               ! Linear system coefficients
668               !
669               zAV(ji,jj) = -   zt11V           * e1v(ji,jj-1) -   zt21V           * r1_e1v(ji,jj-1)
670               zBV(ji,jj) =   ( zt11V + zt12V ) * e1v(ji,jj)   + ( zt21V + zt22V ) * r1_e1v(ji,jj)   + ( zt122V + zt121V ) * r1_e2v(ji,jj)
671               zCV(ji,jj) = -   zt12V           * e1v(ji,jj+1) -   zt22V           * r1_e1v(ji,jj+1)
672         
673               zDV(ji,jj) =     zt121V * r1_e2v(ji-1,jj)
674               zEV(ji,jj) =     zt122V * r1_e2v(ji+1,jj)
675                 
676               !-----------------------------------------------------
677               ! -- Ocean-ice drag and acceleration LHS contribution
678               !-----------------------------------------------------
679               zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj)
680               zBV(ji,jj) = zBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj)
681         
682            END DO
683         END DO
684
685      !------------------------------------------------------------------------------!
686      !
687      ! --- Inner loop: solve linear system, check convergence
688      !
689      !------------------------------------------------------------------------------!
690               
691         ! Inner loop solves the linear problem .. requires 1500 iterations
692         ll_u_iterate = .TRUE.
693         ll_v_iterate = .TRUE.
694
695         DO i_inn = 1, nn_ninn_vp ! inner loop iterations
696
697            !--- mitgcm computes initial value of residual here...
698
699            i_inn_tot             = i_inn_tot + 1
700            ! l_full_nf_update = i_inn_tot == nn_nvp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
701
702            zu_b(:,:)       = u_ice(:,:) ! velocity at previous inner-iterate
703            zv_b(:,:)       = v_ice(:,:)
704
705            IF ( ll_u_iterate .OR. ll_v_iterate )   THEN
706
707                                           ! ---------------------------- !
708               IF ( ll_u_iterate ) THEN    ! --- Solve for u-velocity --- !
709                                           ! ---------------------------- !
710
711                  ! What follows could be subroutinized...
712     
713                  ! Thomas Algorithm  for tridiagonal solver
714                  ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F
715
716                  zFU(:,:)       = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; 
717                 
718                  DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! )
719                 
720                     ! OPT: could be even better optimized with a true red-black SOR
721     
722                     IF ( jn == 1 ) THEN   ;   jj_min = 2 
723                     ELSE                  ;   jj_min = 3
724                     ENDIF
725
726                     DO jj = jj_min, jpj - 1, nn_zebra_vp
727                     
728                        !------------------------
729                        ! Independent term (zFU)
730                        !------------------------
731                        DO ji = 2, jpi - 1   
732                           ! note: these are key lines linking information between processors
733                           ! u_ice/v_ice need to be lbc_linked
734
735                           ! sub-domain boundary condition substitution
736                           ! see Zhang and Hibler, 1997, Appendix B
737                           zAA3 = 0._wp
738                           IF ( ji == 2 )         zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj)
739                           IF ( ji == jpi - 1 )   zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj)
740
741                           ! right hand side
742                           zFU(ji,jj) = ( zrhsu(ji,jj) &                                      ! right-hand side terms
743                               &      +   zAA3         &                                      ! boundary condition translation
744                               &      +   zDU(ji,jj) * u_ice(ji,jj-1)   &                     ! internal force, j-1
745                               &      +   zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1)      ! internal force, j+1
746
747                        END DO
748
749                     END DO
750                     
751                     !---------------
752                     ! Forward sweep
753                     !---------------   
754                     DO jj = jj_min, jpj - 1, nn_zebra_vp
755     
756                        zBU_prime(2,jj)     = zBU(2,jj)
757                        zFU_prime(2,jj)     = zFU(2,jj)
758
759                        DO ji = 3, jpi - 1
760
761                           zfac             = SIGN( 1._wp , zBU(ji-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU(ji-1,jj) ) - epsi20 ) )
762                           zw               = zfac * zAU(ji,jj) / MAX ( ABS( zBU(ji-1,jj) ) , epsi20 ) 
763                           zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj)
764                           zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj)
765
766                        END DO
767
768                     END DO
769                                                                                                     
770                     !-----------------------------
771                     ! Backward sweep & relaxation
772                     !-----------------------------
773
774                     DO jj = jj_min, jpj - 1, nn_zebra_vp
775                   
776                        ! --- Backward sweep
777
778                        ! last row
779                        zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) )
780                        u_ice(jpi-1,jj)    = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 
781                                           &            * umask(jpi-1,jj,1)
782
783                        DO ji = jpi - 2 , 2, -1 ! all other rows    !  ---> original backward loop
784                           zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) )
785                           u_ice(ji,jj)    = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1)   & 
786                                           &                  / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 ) 
787                        END DO
788
789                        !--- Relaxation and masking (for low-ice/no-ice cases)
790                        DO ji = 2, jpi - 1   
791                       
792                           u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) ! relaxation
793                           
794                           u_ice(ji,jj) =   zmsk00x(ji,jj)                                        &   ! masking
795                              &         * (           zmsk01x(ji,jj)   * u_ice(ji,jj)             &
796                              &           + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp     ) * umask(ji,jj,1)
797                           
798                        END DO
799
800                     END DO ! jj
801                     
802                  END DO ! zebra loop
803
804               ENDIF !   ll_u_iterate
805
806               !                           ! ---------------------------- !
807               IF ( ll_v_iterate ) THEN    ! --- Solve for V-velocity --- !
808               !                           ! ---------------------------- !
809                                         
810                  ! MV OPT: what follows could be subroutinized...                 
811                  ! Thomas Algorithm  for tridiagonal solver
812                  ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F
813                  ! It is intentional to have a ji then jj loop for V-velocity
814                  !!! ZH97 explain it is critical for convergence speed
815
816                  zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; 
817
818                  DO jn = 1, nn_zebra_vp ! "zebra" loop
819     
820                     IF ( jn == 1 ) THEN   ;   ji_min = 2 
821                     ELSE                  ;   ji_min = 3
822                     ENDIF
823
824                     DO ji = ji_min, jpi - 1, nn_zebra_vp 
825                     
826                        !------------------------
827                        ! Independent term (zFV)
828                        !------------------------
829                        DO jj = 2, jpj - 1
830
831                           ! subdomain boundary condition substitution (check it is correctly applied !!!)
832                           ! see Zhang and Hibler, 1997, Appendix B
833                           zAA3 = 0._wp
834                           IF ( jj == 2 )       zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1)
835                           IF ( jj == jpj - 1 ) 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) ) * vmask(ji,jj,1)   ! internal force, j+1
842
843                        END DO
844                       
845                     END DO
846
847                     !---------------
848                     ! Forward sweep
849                     !---------------
850                     DO ji = ji_min, jpi - 1, nn_zebra_vp 
851                     
852                        zBV_prime(ji,2)     = zBV(ji,2)
853                        zFV_prime(ji,2)     = zFV(ji,2)
854
855                        DO jj = 3, jpj - 1 
856
857                           zfac             = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) )
858                           zw               = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 )
859                           zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1)
860                           zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 
861
862                        END DO
863
864                     END DO
865
866                     !-----------------------------
867                     ! Backward sweep & relaxation
868                     !-----------------------------
869                     DO ji = ji_min, jpi - 1, nn_zebra_vp 
870                   
871                        ! --- Backward sweep
872                        ! last row
873                        zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) )
874                        v_ice(ji,jpj-1)  = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) & 
875                                         &         * vmask(ji,jpj-1,1)  ! last row
876
877                        ! other rows
878                        DO jj = jpj-2, 2, -1 ! original back loop
879                           zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) )
880                           v_ice(ji,jj)   = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) &
881                                          &  / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 )       
882                        END DO           
883                                                   
884                        ! --- Relaxation & masking
885                        DO jj = 2, jpj - 1
886                       
887                            v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) )    ! relaxation
888                           
889                            v_ice(ji,jj) =   zmsk00y(ji,jj)                                        &      ! masking
890                              &         * (           zmsk01y(ji,jj)   * v_ice(ji,jj)              &
891                              &           + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) * vmask(ji,jj,1)
892
893                        END DO ! jj
894                       
895                     END DO ! ji
896                     
897                  END DO ! zebra loop
898                                   
899               ENDIF !   ll_v_iterate
900
901               ! I suspect the communication should go into the zebra loop if we want reproducibility
902               CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
903                             
904               !--------------------------------------------------------------------------------------
905               ! -- Check convergence based on maximum velocity difference, continue or stop the loop
906               !--------------------------------------------------------------------------------------
907
908               !------
909               ! on U
910               !------
911               ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check
912               ! then we must optimize the use of the mpp_max, which is prohibitive                           
913               zuerr_max  = 0._wp
914                               
915               IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN
916
917                  ! - Maximum U-velocity difference               
918                  zuerr(:,:) = 0._wp
919                  DO jj = 2, jpj - 1
920                     DO ji = 2, jpi - 1
921                        zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 
922                     END DO
923                  END DO
924                  zuerr_max = MAXVAL( zuerr )
925                  CALL mpp_max( 'icedyn_rhg_evp', zuerr_max )   ! max over the global domain - damned!
926
927                  ! - Stop if max :error is too large ("safeguard against bad forcing" of original Zhang routine)
928                  IF ( i_inn > 1 .AND. zuerr_max > rn_uerr_max_vp ) THEN
929                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
930                      ll_u_iterate = .FALSE.
931                  ENDIF
932                 
933                  ! - Stop if error small enough
934                  IF ( zuerr_max < rn_uerr_min_vp ) THEN                                       
935                      IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! "
936                      ll_u_iterate = .FALSE.
937                  ENDIF
938                                               
939               ENDIF ! ll_u_iterate
940
941               !------
942               ! on V
943               !------
944               zverr_max = 0._wp
945               
946               IF ( ll_v_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN
947               
948                  ! - Maximum V-velocity difference
949                  zverr(:,:)   = 0._wp   
950                  DO jj = 2, jpj - 1
951                     DO ji = 2, jpi - 1
952                        zverr(ji,jj) = ABS ( ( v_ice(ji,jj) - zv_b(ji,jj) ) ) * vmask(ji,jj,1)
953                     END DO
954                  END DO
955                  zverr_max = MAXVAL( zverr )
956                  CALL mpp_max( 'icedyn_rhg_evp', zverr_max )   ! max over the global domain - damned!
957
958                  ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)
959                  IF ( i_inn > 1 .AND. zverr_max > rn_uerr_max_vp ) THEN
960                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zverr_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
961                      ll_v_iterate = .FALSE.
962                  ENDIF
963                 
964                  ! - Stop if error small enough
965                  IF ( zverr_max < rn_uerr_min_vp ) THEN                                       
966                      IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! "
967                      ll_v_iterate = .FALSE.
968                  ENDIF
969                 
970               ENDIF ! ll_v_iterate
971
972            ENDIF ! ---    end ll_u_iterate or ll_v_iterate
973               
974            !---------------------------------------------------------------------------------------
975            !
976            ! --- Calculate extra convergence diagnostics and save them
977            !
978            !---------------------------------------------------------------------------------------
979            IF( ln_rhg_chkcvg .AND. MOD ( i_inn - 1, nn_cvgchk_vp ) == 0 ) THEN
980
981               CALL rhg_cvg_vp( kt, i_out, i_inn, i_inn_tot, nn_nout_vp, nn_ninn_vp, nn_nvp,        &
982                      &         u_ice, v_ice, zu_b, zv_b, zu_c, zv_c,                               &
983                      &         zmt, za_iU, za_iV, zuerr_max, zverr_max, zglob_area,                &
984                      &         zrhsu, zAU, zBU, zCU, zDU, zEU, zFU,                                &
985                      &         zrhsv, zAV, zBV, zCV, zDV, zEV, zFV,                                &
986                                zvel_res, zvel_diff )
987
988            ENDIF
989
990         END DO ! i_inn, end of inner loop
991
992      END DO ! End of outer loop (i_out) =============================================================================================
993
994      IF( ln_rhg_chkcvg ) THEN
995         
996         IF( iom_use('velo_res') )   CALL iom_put( 'velo_res', zvel_res  )   ! linear system residual  @last inner&outer iteration
997         IF( iom_use('velo_ero') )   CALL iom_put( 'velo_ero', zvel_diff )   ! abs velocity difference @last outer iteration
998         IF( iom_use('uice_eri') )   CALL iom_put( 'uice_eri', zuerr     )   ! abs velocity difference @last inner iteration
999         IF( iom_use('vice_eri') )   CALL iom_put( 'vice_eri', zverr     )   ! abs velocity difference @last inner iteration
1000
1001         DEALLOCATE( zvel_res , zvel_diff )
1002       
1003      ENDIF ! ln_rhg_chkcvg
1004
1005      !------------------------------------------------------------------------------!
1006      !
1007      ! --- Recompute delta, shear and div (inputs for mechanical redistribution)
1008      !
1009      !------------------------------------------------------------------------------!
1010      !
1011      ! MV OPT: subroutinize ?
1012
1013      DO jj = 1, jpj - 1
1014         DO ji = 1, jpi - 1
1015
1016            ! shear at F points
1017            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
1018               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
1019               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
1020
1021         END DO
1022      END DO           
1023     
1024      DO jj = 2, jpj - 1
1025         DO ji = 2, jpi - 1 !
1026           
1027            ! tension**2 at T points
1028            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
1029               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
1030               &   ) * r1_e1e2t(ji,jj)
1031            zdt2 = zdt * zdt
1032           
1033            ztens(ji,jj)    = zdt
1034           
1035            ! shear**2 at T points (doc eq. A16)
1036            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
1037               &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
1038               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
1039           
1040            ! maximum shear rate at T points (includees tension, output only)
1041            pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) ! i think this is maximum shear rate and not actual shear --- i'm not totally sure here
1042
1043            ! shear at T-points
1044            zshear(ji,jj)   = SQRT( zds2 )
1045
1046            ! divergence at T points
1047            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
1048               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
1049               &             ) * r1_e1e2t(ji,jj)
1050
1051            zdiv2          =  pdivu_i(ji,jj) *  pdivu_i(ji,jj)
1052           
1053            ! delta at T points
1054            zdelta         = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
1055            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
1056
1057            pdelta_i(ji,jj) = zdelta + rn_creepl ! * rswitch
1058
1059         END DO
1060      END DO
1061     
1062      CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
1063     
1064      !------------------------------------------------------------------------------!
1065      !
1066      ! --- Diagnostics
1067      !
1068      !------------------------------------------------------------------------------!
1069      !
1070      ! MV OPT: subroutinize ?
1071      !
1072      !----------------------------------
1073      ! --- Recompute stresses if needed
1074      !----------------------------------
1075      !
1076      ! ---- Sea ice stresses at T-points
1077      IF ( iom_use('normstr')    .OR. iom_use('sheastr')    .OR. &
1078     &     iom_use('intstrx')    .OR. iom_use('intstry')    .OR. &
1079     &     iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
1080     
1081         DO jj = 2, jpj - 1
1082            DO ji = 2, jpi - 1
1083               zp_delstar_t(ji,jj)     =   strength(ji,jj) / pdelta_i(ji,jj) 
1084               zfac                    =   zp_delstar_t(ji,jj) 
1085               zs1(ji,jj)              =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
1086               zs2(ji,jj)              =   zfac * z1_ecc2 * ztens(ji,jj)
1087               zs12(ji,jj)             =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bug 12 nov
1088            END DO
1089         END DO
1090
1091         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. )
1092     
1093      ENDIF
1094     
1095      ! ---- s12 at F-points     
1096      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1097
1098         DO jj = 1, jpj - 1
1099            DO ji = 1, jpi - 1
1100           
1101               ! P/delta* at F points
1102               zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) )
1103               
1104               ! s12 at F-points
1105               zs12f(ji,jj) = zp_delstar_f * z1_ecc2 * zds(ji,jj)
1106               
1107            END DO
1108         END DO
1109
1110         CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. )
1111         
1112      ENDIF
1113
1114      !
1115      !-----------------------
1116      ! --- Store diagnostics
1117      !-----------------------
1118      !
1119      ! --- Ice-ocean, ice-atm. & ice-ocean bottom (landfast) stresses --- !
1120      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
1121         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
1122
1123         ALLOCATE( ztaux_oi(jpi,jpj) , ztauy_oi(jpi,jpj) )
1124
1125         !--- Recalculate oceanic stress at last inner iteration
1126         DO jj = 2, jpj - 1
1127            DO ji = 2, jpi - 1
1128
1129                !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
1130                zu_cV            = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)
1131                zv_cU            = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)
1132               
1133                !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
1134                zCwU(ji,jj)          = za_iU(ji,jj) * zrhoco * SQRT( ( u_ice(ji,jj) - u_oce (ji,jj) ) * ( u_ice(ji,jj) - u_oce (ji,jj) )  &
1135                  &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
1136                zCwV(ji,jj)          = za_iV(ji,jj) * zrhoco * SQRT( ( v_ice(ji,jj) - v_oce (ji,jj) ) * ( v_ice(ji,jj) - v_oce (ji,jj) )  &
1137                  &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
1138                 
1139                !--- Ocean-ice stress
1140                ztaux_oi(ji,jj) = zCwU(ji,jj) * ( u_oce(ji,jj) - u_ice(ji,jj) )
1141                ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) )
1142               
1143            END DO
1144         END DO
1145         
1146         !
1147         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, &
1148!            &                                 ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
1149         !
1150         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
1151         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
1152         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
1153         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
1154!        CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
1155!        CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
1156
1157         DEALLOCATE( ztaux_oi , ztauy_oi )
1158
1159      ENDIF
1160       
1161      ! --- Divergence, shear and strength --- !
1162      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
1163      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! maximum shear rate
1164      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
1165      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
1166
1167      ! --- Stress tensor invariants (SIMIP diags) --- !
1168      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
1169         !
1170         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags.
1171         ! Definitions following Coon (1974) and Feltham (2008)
1172         !
1173         ! sigma1, sigma2, sigma12 are useful (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
1174         ! however these are NOT stress tensor components, neither stress invariants, nor stress principal components
1175         !
1176         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
1177         !         
1178         DO jj = 2, jpj - 1
1179            DO ji = 2, jpi - 1
1180               ! Stress invariants
1181               zsig_I(ji,jj)    =   zs1(ji,jj) * 0.5_wp   ! 1st invariant, aka average normal stress aka negative pressure
1182               zsig_II(ji,jj)   =   0.5_wp * SQRT ( zs2(ji,jj) * zs2(ji,jj) + 4. * zs12(ji,jj) * zs12(ji,jj) )   ! 2nd invariant, aka maximum shear stress               
1183            END DO
1184         END DO
1185
1186         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
1187         
1188         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,   zsig_I(:,:)  * zmsk00(:,:) ) ! Normal stress
1189         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' ,   zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
1190         
1191         DEALLOCATE ( zsig_I, zsig_II )
1192         
1193      ENDIF
1194
1195      ! --- Normalized stress tensor principal components --- !
1196      ! These are used to plot the normalized yield curve (Lemieux & Dupont, GMD 2020)
1197      ! To plot the yield curve and evaluate physical convergence, they have two recommendations
1198      ! Recommendation 1 : Use ice strength, not replacement pressure
1199      ! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765)
1200      ! R2 means we need to recompute stresses
1201
1202      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
1203         !
1204         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
1205         !         
1206         DO jj = 2, jpj - 1
1207            DO ji = 2, jpi - 1
1208               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
1209               !                        and **deformations** at current iterates
1210               !                        following Lemieux & Dupont (2020)
1211               zfac             =   zp_delstar_t(ji,jj)
1212               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltat(ji,jj) )
1213               zsig2            =   zfac * z1_ecc2 * ztens(ji,jj)
1214               zsig12           =   zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp ! Bugfix 12 Nov
1215               
1216               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
1217               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                       ! 1st invariant
1218               zsig_II(ji,jj)   =   0.5_wp * SQRT ( zsig2 * zsig2 + 4. *zsig12 * zsig12 )   ! 2nd invariant
1219
1220               ! Normalized  principal stresses (used to display the ellipse)
1221               z1_strength      =   1._wp / MAX ( 1._wp , strength(ji,jj) )
1222               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
1223               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
1224            END DO
1225         END DO
1226         !
1227         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
1228         !
1229         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
1230         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
1231
1232         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II )
1233         
1234      ENDIF
1235
1236      ! --- SIMIP, terms of tendency for momentum equation  --- !
1237      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
1238         & iom_use('corstrx') .OR. iom_use('corstry') ) THEN
1239
1240         ! --- Recalculate Coriolis stress at last inner iteration
1241         DO jj = 2, jpj - 1
1242            DO ji = 2, jpi - 1
1243                ! --- U-component
1244                zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  &
1245                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
1246                           &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
1247                zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  &
1248                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
1249                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
1250            END DO
1251         END DO
1252         !
1253         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., &
1254            &                                 zCorU, 'U', -1., zCorV, 'V', -1. )
1255         !
1256         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
1257         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
1258         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
1259         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
1260
1261      ENDIF
1262     
1263      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1264
1265         ! Recalculate internal forces (divergence of stress tensor) at last inner iteration
1266         DO jj = 2, jpj - 1
1267            DO ji = 2, jpi - 1
1268
1269               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
1270                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
1271                  &                    ) * r1_e2u(ji,jj)                                                                      &
1272                  &                  + ( zs12f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
1273                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
1274                  &                  ) * r1_e1e2u(ji,jj)
1275
1276               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
1277                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
1278                  &                    ) * r1_e1v(ji,jj)                                                                      &
1279                  &                  + ( zs12f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
1280                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
1281                  &                  ) * r1_e1e2v(ji,jj)
1282
1283            END DO
1284         END DO
1285           
1286         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. )
1287         
1288         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
1289         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
1290         
1291      ENDIF
1292
1293      ! --- Ice & snow mass and ice area transports
1294      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
1295         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
1296         !
1297         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
1298            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
1299         !
1300         DO jj = 2, jpj - 1
1301            DO ji = 2, jpi - 1
1302               ! 2D ice mass, snow mass, area transport arrays (X, Y)
1303               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
1304               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
1305
1306               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
1307               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
1308
1309               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
1310               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
1311
1312               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
1313               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
1314
1315            END DO
1316         END DO
1317
1318         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
1319            &                                 zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
1320            &                                 zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
1321
1322         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
1323         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport
1324         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
1325         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
1326         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
1327         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
1328
1329         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
1330            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
1331
1332      ENDIF
1333
1334      DEALLOCATE( zmsk00, zmsk15 )
1335
1336   END SUBROUTINE ice_dyn_rhg_vp
1337   
1338   
1339   SUBROUTINE rhg_cvg_vp( kt, kitout, kitinn, kitinntot, kitoutmax, kitinnmax, kitinntotmax , &
1340                  &       pu, pv, pub, pvb, pub_outer, pvb_outer                     , &
1341                  &       pmt, pat_iu, pat_iv, puerr_max, pverr_max, pglob_area      , &
1342                  &       prhsu, pAU, pBU, pCU, pDU, pEU, pFU                        , &
1343                  &       prhsv, pAV, pBV, pCV, pDV, pEV, pFV                        , &   
1344                  &       pvel_res, pvel_diff                                            )
1345      !!
1346      !!----------------------------------------------------------------------
1347      !!                    ***  ROUTINE rhg_cvg_vp  ***
1348      !!                     
1349      !! ** Purpose :   check convergence of VP ice rheology
1350      !!
1351      !! ** Method  :   create a file ice_cvg.nc containing a few convergence diagnostics
1352      !!                This routine is called every sub-iteration, so it is cpu expensive
1353      !!
1354      !!                Calculates / stores
1355      !!                   - maximum absolute U-V difference (uice_cvg, u_dif, v_dif, m/s)
1356      !!                   - 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)
1357      !!                   - mean kinetic energy (mke_ice, J/m2)
1358      !!
1359      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)   
1360      !!
1361      !!----------------------------------------------------------------------
1362      !!
1363      INTEGER ,                 INTENT(in) ::   kt, kitout, kitinn, kitinntot    ! ocean model iterate, outer, inner and total n-iterations
1364      INTEGER ,                 INTENT(in) ::   kitoutmax, kitinnmax             ! max number of outer & inner iterations
1365      INTEGER ,                 INTENT(in) ::   kitinntotmax                     ! max number of total sub-iterations
1366      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb                 ! now & sub-iter-before velocities
1367      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pub_outer, pvb_outer             ! velocities @before outer iterations
1368      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pmt, pat_iu, pat_iv              ! mass at T-point, ice concentration at U&V
1369      REAL(wp),                 INTENT(in) ::   puerr_max, pverr_max             ! absolute mean velocity difference
1370      REAL(wp),                 INTENT(in) ::   pglob_area                       ! global ice area
1371      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsu, pAU, pBU, pCU, pDU, pEU, pFU ! linear system coefficients
1372      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsv, pAV, pBV, pCV, pDV, pEV, pFV
1373      REAL(wp), DIMENSION(:,:), INTENT(inout) ::  pvel_res                       ! velocity residual @last inner iteration
1374      REAL(wp), DIMENSION(:,:), INTENT(inout) ::  pvel_diff                      ! velocity difference @last outer iteration
1375      !!
1376
1377      INTEGER           ::   idtime, istatus, ix_dim, iy_dim
1378      INTEGER           ::   ji, jj          ! dummy loop indices
1379      INTEGER           ::   it_inn_file, it_out_file
1380      REAL(wp)          ::   zu_res_mean, zv_res_mean, zvel_res_mean                  ! mean residuals of the linear system
1381      REAL(wp)          ::   zu_mad, zv_mad, zvel_mad                                 ! mean absolute deviation, sub-iterates
1382      REAL(wp)          ::   zu_mad_outer, zv_mad_outer, zvel_mad_outer               ! mean absolute deviation, outer-iterates
1383      REAL(wp)          ::   zvel_err_max, zmke, zu, zv                               ! local scalars
1384      REAL(wp)          ::   z1_pglob_area                                            ! inverse global ice area
1385
1386      REAL(wp), DIMENSION(jpi,jpj) ::   zu_res, zv_res, zvel2                         ! local arrays
1387      REAL(wp), DIMENSION(jpi,jpj) ::   zu_diff, zv_diff                              ! local arrays
1388                                                                             
1389      CHARACTER(len=20) ::   clname
1390      !!----------------------------------------------------------------------
1391
1392
1393      IF( lwp ) THEN
1394
1395         WRITE(numout,*)
1396         WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control'
1397         WRITE(numout,*) '~~~~~~~~~~~'
1398         WRITE(numout,*) ' kt          =  : ', kt
1399         WRITE(numout,*) ' kitout      =  : ', kitout
1400         WRITE(numout,*) ' kitinn      =  : ', kitinn
1401         WRITE(numout,*) ' kitinntot   =  : ', kitinntot
1402         WRITE(numout,*) ' kitoutmax (nn_nout_vp) =  ', kitoutmax
1403         WRITE(numout,*) ' kitinnmax (nn_ninn_vp) =  ', kitinnmax
1404         WRITE(numout,*) ' kitinntotmax (nn_nvp)  =  ', kitinntotmax
1405         WRITE(numout,*)
1406
1407      ENDIF
1408
1409      z1_pglob_area = 1._wp / pglob_area      ! inverse global ice area
1410
1411      ! create file
1412      IF( kt == nit000 .AND. kitinntot == 1 ) THEN
1413         !
1414         IF( lwm ) THEN
1415
1416            clname = 'ice_cvg.nc'
1417            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
1418            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
1419
1420            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
1421            istatus = NF90_DEF_DIM( ncvgid, 'x'     , jpi, ix_dim )
1422            istatus = NF90_DEF_DIM( ncvgid, 'y'     , jpj, iy_dim )
1423
1424            istatus = NF90_DEF_VAR( ncvgid, 'u_res'         , NF90_DOUBLE  , (/ idtime /), nvarid_ures )
1425            istatus = NF90_DEF_VAR( ncvgid, 'v_res'         , NF90_DOUBLE  , (/ idtime /), nvarid_vres )
1426            istatus = NF90_DEF_VAR( ncvgid, 'vel_res'       , NF90_DOUBLE  , (/ idtime /), nvarid_velres )
1427
1428            istatus = NF90_DEF_VAR( ncvgid, 'uerr_max_sub'  , NF90_DOUBLE  , (/ idtime /), nvarid_uerr_max )
1429            istatus = NF90_DEF_VAR( ncvgid, 'verr_max_sub'  , NF90_DOUBLE  , (/ idtime /), nvarid_verr_max )
1430            istatus = NF90_DEF_VAR( ncvgid, 'velerr_max_sub', NF90_DOUBLE  , (/ idtime /), nvarid_velerr_max )
1431
1432            istatus = NF90_DEF_VAR( ncvgid, 'umad_sub'      , NF90_DOUBLE  , (/ idtime /), nvarid_umad )
1433            istatus = NF90_DEF_VAR( ncvgid, 'vmad_sub'      , NF90_DOUBLE  , (/ idtime /), nvarid_vmad )
1434            istatus = NF90_DEF_VAR( ncvgid, 'velmad_sub'    , NF90_DOUBLE  , (/ idtime /), nvarid_velmad )
1435           
1436            istatus = NF90_DEF_VAR( ncvgid, 'umad_outer'    , NF90_DOUBLE  , (/ idtime /), nvarid_umad_outer   )
1437            istatus = NF90_DEF_VAR( ncvgid, 'vmad_outer'    , NF90_DOUBLE  , (/ idtime /), nvarid_vmad_outer   )
1438            istatus = NF90_DEF_VAR( ncvgid, 'velmad_outer'  , NF90_DOUBLE  , (/ idtime /), nvarid_velmad_outer )
1439
1440            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE  , (/ idtime /), nvarid_mke )
1441
1442            istatus = NF90_ENDDEF(ncvgid)
1443
1444         ENDIF
1445         !
1446      ENDIF
1447
1448      !------------------------------------------------------------
1449      !
1450      ! Max absolute velocity difference with previous sub-iterate
1451      ! ( zvel_err_max )
1452      !
1453      !------------------------------------------------------------
1454      !
1455      ! This comes from the criterion used to stop the iterative procedure
1456      zvel_err_max   = 0.5_wp * ( puerr_max + pverr_max ) ! average of U- and V- maximum error over the whole domain
1457
1458      !----------------------------------------------
1459      !
1460      ! Mean-absolute-deviation (sub-iterates)
1461      ! ( zu_mad, zv_mad, zvel_mad)
1462      !
1463      !----------------------------------------------
1464      !
1465      ! U
1466      zu_diff(:,:) = 0._wp
1467      DO jj = 2, jpj - 1
1468         DO ji = 2, jpi - 1
1469            zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * z1_pglob_area
1470         END DO
1471      END DO
1472
1473      ! V
1474      zv_diff(:,:)   = 0._wp
1475      DO jj = 2, jpj - 1
1476         DO ji = 2, jpi - 1
1477            zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * z1_pglob_area
1478         END DO
1479      END DO
1480
1481
1482      ! global sum & U-V average
1483      CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff,  'U',  1., zv_diff , 'V',  1. )
1484      zu_mad   = glob_sum( 'icedyn_rhg_vp : ', zu_diff )
1485      zv_mad   = glob_sum( 'icedyn_rhg_vp : ', zv_diff )
1486
1487      zvel_mad = 0.5_wp * ( zu_mad + zv_mad )
1488
1489      !-----------------------------------------------
1490      !
1491      ! Mean-absolute-deviation (outer-iterates)
1492      ! ( zu_mad_outer, zv_mad_outer, zvel_mad_outer)
1493      !
1494      !-----------------------------------------------
1495      !
1496      IF ( kitinn == kitinnmax ) THEN ! only work at the end of outer iterates
1497
1498         ! * U
1499         zu_diff(:,:) = 0._wp
1500         DO jj = 2, jpj - 1
1501            DO ji = 2, jpi - 1
1502               zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub_outer(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * &
1503                              &    z1_pglob_area
1504            END DO
1505         END DO
1506
1507         ! * V
1508         zv_diff(:,:)   = 0._wp
1509         DO jj = 2, jpj - 1
1510            DO ji = 2, jpi - 1
1511               zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb_outer(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * &
1512                              &    z1_pglob_area
1513            END DO
1514         END DO
1515
1516         ! Global ice-concentration, grid-cell-area weighted mean
1517         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff,  'U',  1., zv_diff , 'V',  1. ) ! abs behaves as a scalar no ?
1518
1519         zu_mad_outer   = glob_sum( 'icedyn_rhg_vp : ', zu_diff )
1520         zv_mad_outer   = glob_sum( 'icedyn_rhg_vp : ', zv_diff )
1521   
1522         ! Average of both U & V
1523         zvel_mad_outer = 0.5_wp * ( zu_mad_outer + zv_mad_outer )
1524                 
1525      ENDIF
1526
1527      ! --- Spatially-resolved absolute difference to send back to main routine
1528      ! (last iteration only, T-point)
1529
1530      IF ( kitinntot == kitinntotmax) THEN
1531
1532         zu_diff(:,:) = 0._wp
1533         zv_diff(:,:) = 0._wp
1534
1535         DO jj = 2, jpj - 1
1536            DO ji = 2, jpi - 1
1537
1538               zu_diff(ji,jj) = ( ABS ( ( pu(ji-1,jj) - pub_outer(ji-1,jj) ) ) * umask(ji-1,jj,1) &
1539    &                           + ABS ( ( pu(ji,jj  ) - pub_outer(ji,jj)   ) ) * umask(ji,jj,1) ) &
1540    &                           / ( umask(ji-1,jj,1) + umask(ji,jj,1) )
1541
1542               zv_diff(ji,jj) = ( ABS ( ( pv(ji,jj-1)   - pvb_outer(ji,jj-1) ) ) * vmask(ji,jj-1,1) &
1543    &                           + ABS ( ( pv(ji,jj  ) - pvb_outer(ji,jj)     ) ) * vmask(ji,jj,1)   &
1544    &                           / ( vmask(ji,jj-1,1) + vmask(ji,jj,1) ) )
1545     
1546            END DO
1547         END DO
1548
1549         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff,  'T',  1., zv_diff , 'T',  1. )
1550         pvel_diff(:,:) = 0.5_wp * ( zu_diff(:,:) + zv_diff(:,:) )
1551
1552      ELSE
1553
1554         pvel_diff(:,:) = 0._wp
1555
1556      ENDIF
1557
1558      !---------------------------------------
1559      !
1560      ! ---  Mean residual & kinetic energy
1561      !
1562      !---------------------------------------
1563
1564      IF ( kitinntot == 1 ) THEN
1565
1566         zu_res_mean   = 0._wp
1567         zv_res_mean   = 0._wp
1568         zvel_res_mean = 0._wp
1569         zmke          = 0._wp
1570
1571      ELSE
1572
1573         ! * Mean residual (N/m2)
1574         ! Here we take the residual of the linear system (N/m2),
1575         ! We define it as in mitgcm: global area-weighted mean of square-root residual
1576         ! Local residual r = Ax - B expresses to which extent the momentum balance is verified
1577         ! i.e., how close we are to a solution
1578         zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp
1579
1580         DO jj = 2, jpj - 1
1581            DO ji = 2, jpi - 1                                     
1582               zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               &
1583                  &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
1584               zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               &
1585                  &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
1586
1587!              zu_res(ji,jj)  = pFU(ji,jj) - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj)
1588!              zv_res(ji,jj)  = pFV(ji,jj) - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1)
1589   
1590               zu_res(ji,jj)  = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * pat_iu(ji,jj) * e1e2u(ji,jj) * z1_pglob_area
1591               zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * pat_iv(ji,jj) * e1e2v(ji,jj) * z1_pglob_area
1592   
1593            END DO
1594         END DO                 
1595
1596         ! Global ice-concentration, grid-cell-area weighted mean
1597         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_res,  'U', 1., zv_res , 'V', 1. )
1598   
1599         zu_res_mean   = glob_sum( 'ice_rhg_vp', zu_res(:,:) )
1600         zv_res_mean   = glob_sum( 'ice_rhg_vp', zv_res(:,:) )
1601         zvel_res_mean = 0.5_wp * ( zu_res_mean + zv_res_mean )
1602   
1603         ! --- Global mean kinetic energy per unit area (J/m2)
1604         zvel2(:,:) = 0._wp
1605         DO jj = 2, jpj - 1
1606            DO ji = 2, jpi - 1                   
1607               zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point
1608               zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) )
1609               zvel2(ji,jj)  = zu*zu + zv*zv              ! square of ice velocity at T-point 
1610            END DO
1611         END DO
1612         
1613         zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area
1614   
1615      ENDIF ! kitinntot
1616
1617      !--- Spatially-resolved residual at last iteration to send back to main routine (last iteration only)
1618      !--- Calculation @T-point
1619
1620      IF ( kitinntot == kitinntotmax) THEN
1621
1622         DO jj = 2, jpj - 1 ! @ residuals at u and v points
1623            DO ji = 2, jpi - 1
1624
1625               zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               &
1626                  &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
1627               zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               &
1628                  &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
1629
1630               zu_res(ji,jj)  = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) 
1631               zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) 
1632
1633            END DO
1634         END DO
1635
1636         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_res,  'U',  1., zv_res , 'V',  1. )
1637
1638         DO jj = 2, jpj - 1         ! average @T-point
1639            DO ji = 2, jpi - 1
1640               pvel_res(ji,jj) = 0.25_wp * ( zu_res(ji-1,jj) + zu_res(ji,jj) + zv_res(ji,jj-1) + zv_res(ji,jj) )
1641            END DO
1642         END DO
1643
1644         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. )
1645
1646      ELSE
1647
1648         pvel_res(:,:) = 0._wp
1649
1650      ENDIF
1651                 
1652      !                                                ! ==================== !
1653
1654      it_inn_file =  ( kt - nit000 ) * kitinntotmax + kitinntot ! time step in the file
1655      it_out_file =  ( kt - nit000 ) * kitoutmax    + kitout
1656
1657      ! write variables
1658      IF( lwm ) THEN
1659
1660         istatus = NF90_PUT_VAR( ncvgid, nvarid_ures  , (/zu_res_mean/), (/it_inn_file/), (/1/) )        ! Residuals of the linear system, area weighted mean
1661         istatus = NF90_PUT_VAR( ncvgid, nvarid_vres  , (/zv_res_mean/), (/it_inn_file/), (/1/) )        !
1662         istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvel_res_mean/), (/it_inn_file/), (/1/) )      !
1663
1664         istatus = NF90_PUT_VAR( ncvgid, nvarid_uerr_max  , (/puerr_max/), (/it_inn_file/), (/1/) )      ! Max velocit_inn_filey error, sub-it_inn_fileerates
1665         istatus = NF90_PUT_VAR( ncvgid, nvarid_verr_max  , (/pverr_max/), (/it_inn_file/), (/1/) )      !
1666         istatus = NF90_PUT_VAR( ncvgid, nvarid_velerr_max, (/zvel_err_max/), (/it_inn_file/), (/1/) )   !
1667
1668         istatus = NF90_PUT_VAR( ncvgid, nvarid_umad    , (/zu_mad/)  , (/it_inn_file/), (/1/) )         ! velocit_inn_filey MAD, area/sic-weighted, sub-it_inn_fileerates
1669         istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad    , (/zv_mad/)  , (/it_inn_file/), (/1/) )         !
1670         istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad  , (/zvel_mad/), (/it_inn_file/), (/1/) )         !
1671
1672         istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/kitinntot/), (/1/) )                    ! mean kinetic energy
1673
1674         IF ( kitinn == kitinnmax ) THEN ! only print outer mad at the end of inner loop
1675
1676            istatus = NF90_PUT_VAR( ncvgid, nvarid_umad_outer    , (/zu_mad_outer/)  , (/it_out_file/), (/1/) )   ! velocity MAD, area/sic-weighted, outer-iterates
1677            istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad_outer    , (/zv_mad_outer/)  , (/it_out_file/), (/1/) )   !
1678            istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad_outer  , (/zvel_mad_outer/), (/it_out_file/), (/1/) )   !
1679
1680         ENDIF
1681
1682         IF( kt == nitend - nn_fsbc + 1 .AND. kitinntot == kitinntotmax )    istatus = NF90_CLOSE( ncvgid )
1683      ENDIF
1684     
1685   END SUBROUTINE rhg_cvg_vp
1686   
1687
1688   
1689#else
1690   !!----------------------------------------------------------------------
1691   !!   Default option         Empty module           NO SI3 sea-ice model
1692   !!----------------------------------------------------------------------
1693#endif
1694
1695   !!==============================================================================
1696END MODULE icedyn_rhg_vp
1697
Note: See TracBrowser for help on using the repository browser.