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 @ 15515

Last change on this file since 15515 was 15515, checked in by vancop, 8 months ago

some bugfixes on diagnostics

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