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

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

northfold and stress diagnostic issues in VP rheology

File size: 97.4 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) ::   zten_i                          ! Tension
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)
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            zten_i(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            ! shear at T points
1229            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
1230
1231            ! divergence at T points
1232            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
1233               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
1234               &             ) * r1_e1e2t(ji,jj)
1235           
1236            ! delta at T points
1237            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
1238            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
1239           
1240            !pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
1241            pdelta_i(ji,jj) = zdelta + rn_creepl
1242
1243         END DO
1244      END DO
1245
1246      IF ( lwp ) WRITE(numout,*) ' Deformation recalculated '
1247     
1248      CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
1249     
1250      !------------------------------------------------------------------------------!
1251      !
1252      ! --- Diagnostics
1253      !
1254      !------------------------------------------------------------------------------!
1255      !
1256      ! MV OPT: subroutinize ?
1257      !
1258      !----------------------------------
1259      ! --- Recompute stresses if needed
1260      !----------------------------------
1261      !
1262      ! ---- Sea ice stresses at T-points
1263      IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1264     
1265         DO jj = 2, jpj - 1
1266            DO ji = 2, jpi - 1
1267               zp_deltastar_t(ji,jj)   =   strength(ji,jj) / pdelta_i(ji,jj) 
1268               zfac                    =   zp_deltastar_t(ji,jj) 
1269               zs1(ji,jj)              =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
1270               zs2(ji,jj)              =   zfac * z1_ecc2 * zten_i(ji,jj)
1271               zs12(ji,jj)             =   zfac * z1_ecc2 * pshear_i(ji,jj)
1272            END DO
1273         END DO
1274
1275         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. )
1276     
1277      ENDIF
1278     
1279      ! ---- s12 at F-points     
1280      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1281
1282         DO jj = 1, jpj - 1
1283            DO ji = 1, jpi - 1
1284           
1285               ! P/delta* at F points
1286               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) )
1287               
1288               ! s12 at F-points
1289               zs12f(ji,jj) = zp_deltastar_f * z1_ecc2 * zds(ji,jj)
1290               
1291            END DO
1292         END DO
1293
1294         CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. )
1295         
1296      ENDIF
1297
1298      IF ( lwp ) WRITE(numout,*) ' zs12f recalculated '
1299
1300      !
1301      !-----------------------
1302      ! --- Store diagnostics
1303      !-----------------------
1304      !
1305      ! --- Ice-ocean, ice-atm. & ice-ocean bottom (landfast) stresses --- !
1306      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
1307         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
1308
1309         ALLOCATE( ztaux_oi(jpi,jpj) , ztauy_oi(jpi,jpj) )
1310
1311         !--- Recalculate oceanic stress at last inner iteration
1312         DO jj = 2, jpj - 1
1313            DO ji = 2, jpi - 1
1314
1315                !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
1316                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)
1317                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)
1318               
1319                !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
1320                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) )  &
1321                  &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
1322                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) )  &
1323                  &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
1324                 
1325                !--- Ocean-ice stress
1326                ztaux_oi(ji,jj) = zCwU(ji,jj) * ( u_oce(ji,jj) - u_ice(ji,jj) )
1327                ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) )
1328               
1329            END DO
1330         END DO
1331         
1332         !
1333         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, &
1334!            &                                 ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
1335         !
1336         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
1337         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
1338         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
1339         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
1340!        CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
1341!        CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
1342
1343         DEALLOCATE( ztaux_oi , ztauy_oi )
1344
1345      ENDIF
1346       
1347      ! --- Divergence, shear and strength --- !
1348      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
1349      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
1350      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
1351      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
1352
1353      IF ( lwp ) WRITE(numout,*) 'Some terms recalculated '
1354
1355      ! --- Stress tensor invariants (SIMIP diags) --- !
1356      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
1357         !
1358         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags.
1359         ! Definitions following Coon (1974) and Feltham (2008)
1360         !
1361         ! sigma1, sigma2, sigma12 are useful (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
1362         ! however these are NOT stress tensor components, neither stress invariants, nor stress principal components
1363         !
1364         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
1365         !         
1366         DO jj = 2, jpj - 1
1367            DO ji = 2, jpi - 1
1368               ! Stress invariants
1369               zsig_I(ji,jj)    =   zs1(ji,jj) * 0.5_wp                                                   ! 1st invariant, aka average normal stress aka negative pressure
1370               zsig_II(ji,jj)   =   SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj)*zs12(ji,jj) )  ! 2nd invariant, aka maximum shear stress               
1371            END DO
1372         END DO
1373
1374         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
1375         
1376         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,   zsig_I(:,:)  * zmsk00(:,:) ) ! Normal stress
1377         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' ,   zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
1378         
1379         DEALLOCATE ( zsig_I, zsig_II )
1380         
1381      ENDIF
1382
1383      IF ( lwp ) WRITE(numout,*) 'SIMIP work done'
1384
1385      ! --- Normalized stress tensor principal components --- !
1386      ! These are used to plot the normalized yield curve (Lemieux & Dupont, GMD 2020)
1387      ! To plot the yield curve and evaluate physical convergence, they have two recommendations
1388      ! Recommendation 1 : Use ice strength, not replacement pressure
1389      ! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765)
1390      ! R2 means we need to recompute stresses
1391
1392      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
1393         !
1394         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
1395         !         
1396         DO jj = 2, jpj - 1
1397            DO ji = 2, jpi - 1
1398               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
1399               !                        and **deformations** at current iterates
1400               !                        following Lemieux & Dupont (2020)
1401               zfac             =   zp_deltastar_t(ji,jj)
1402               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) )
1403               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
1404               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
1405               
1406               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
1407               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                              ! 1st invariant
1408               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 )   ! 2nd invariant
1409
1410               ! Normalized  principal stresses (used to display the ellipse)
1411               z1_strength      =   1._wp / MAX ( 1._wp , strength(ji,jj) )
1412               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
1413               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
1414            END DO
1415         END DO
1416         IF ( lwp ) WRITE(numout,*) 'Some shitty stress work done'
1417         !
1418         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
1419         !     
1420         IF ( lwp ) WRITE(numout,*) ' Beauaaaarflblbllll '
1421         !
1422         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
1423         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
1424
1425         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II )
1426
1427         IF ( lwp ) WRITE(numout,*) ' So what ??? '
1428         
1429      ENDIF
1430
1431      ! --- SIMIP, terms of tendency for momentum equation  --- !
1432      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
1433         & iom_use('corstrx') .OR. iom_use('corstry') ) THEN
1434
1435         ! --- Recalculate Coriolis stress at last inner iteration
1436         DO jj = 2, jpj - 1
1437            DO ji = 2, jpi - 1
1438                ! --- U-component
1439                zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  &
1440                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
1441                           &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
1442                zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  &
1443                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
1444                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
1445            END DO
1446         END DO
1447         !
1448         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., &
1449            &                                 zCorU, 'U', -1., zCorV, 'V', -1. )
1450         !
1451         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
1452         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
1453         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
1454         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
1455
1456      ENDIF
1457     
1458      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1459
1460         ! Recalculate internal forces (divergence of stress tensor) at last inner iteration
1461         DO jj = 2, jpj - 1
1462            DO ji = 2, jpi - 1
1463
1464               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
1465                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
1466                  &                    ) * r1_e2u(ji,jj)                                                                      &
1467                  &                  + ( zs12f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
1468                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
1469                  &                  ) * r1_e1e2u(ji,jj)
1470
1471               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
1472                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
1473                  &                    ) * r1_e1v(ji,jj)                                                                      &
1474                  &                  + ( zs12f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
1475                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
1476                  &                  ) * r1_e1e2v(ji,jj)
1477
1478            END DO
1479         END DO
1480           
1481         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. )
1482         
1483         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
1484         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
1485         
1486      ENDIF
1487
1488      ! --- Ice & snow mass and ice area transports
1489      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
1490         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
1491         !
1492         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
1493            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
1494         !
1495         DO jj = 2, jpj - 1
1496            DO ji = 2, jpi - 1
1497               ! 2D ice mass, snow mass, area transport arrays (X, Y)
1498               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
1499               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
1500
1501               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
1502               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
1503
1504               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
1505               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
1506
1507               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
1508               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
1509
1510            END DO
1511         END DO
1512
1513         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
1514            &                                 zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
1515            &                                 zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
1516
1517         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
1518         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport
1519         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
1520         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
1521         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
1522         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
1523
1524         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
1525            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
1526
1527      ENDIF
1528
1529      DEALLOCATE( zmsk00, zmsk15 )
1530
1531   END SUBROUTINE ice_dyn_rhg_vp
1532   
1533   
1534   SUBROUTINE rhg_cvg_vp( kt, kitout, kitinn, kitinntot, kitoutmax, kitinnmax, kitinntotmax , &
1535                  &       pu, pv, pub, pvb, pub_outer, pvb_outer                     , &
1536                  &       pmt, pat_iu, pat_iv, puerr_max, pverr_max, pglob_area      , &
1537                  &       prhsu, pAU, pBU, pCU, pDU, pEU, pFU                        , &
1538                  &       prhsv, pAV, pBV, pCV, pDV, pEV, pFV                        , &   
1539                  &       pvel_res, pvel_diff                                            )
1540      !!
1541      !!----------------------------------------------------------------------
1542      !!                    ***  ROUTINE rhg_cvg_vp  ***
1543      !!                     
1544      !! ** Purpose :   check convergence of VP ice rheology
1545      !!
1546      !! ** Method  :   create a file ice_cvg.nc containing a few convergence diagnostics
1547      !!                This routine is called every sub-iteration, so it is cpu expensive
1548      !!
1549      !!                Calculates / stores
1550      !!                   - maximum absolute U-V difference (uice_cvg, u_dif, v_dif, m/s)
1551      !!                   - 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)
1552      !!                   - mean kinetic energy (mke_ice, J/m2)
1553      !!
1554      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)   
1555      !!
1556      !!----------------------------------------------------------------------
1557      !!
1558      INTEGER ,                 INTENT(in) ::   kt, kitout, kitinn, kitinntot    ! ocean model iterate, outer, inner and total n-iterations
1559      INTEGER ,                 INTENT(in) ::   kitoutmax, kitinnmax             ! max number of outer & inner iterations
1560      INTEGER ,                 INTENT(in) ::   kitinntotmax                     ! max number of total sub-iterations
1561      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb                 ! now & sub-iter-before velocities
1562      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pub_outer, pvb_outer             ! velocities @before outer iterations
1563      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pmt, pat_iu, pat_iv              ! mass at T-point, ice concentration at U&V
1564      REAL(wp),                 INTENT(in) ::   puerr_max, pverr_max             ! absolute mean velocity difference
1565      REAL(wp),                 INTENT(in) ::   pglob_area                       ! global ice area
1566      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsu, pAU, pBU, pCU, pDU, pEU, pFU ! linear system coefficients
1567      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsv, pAV, pBV, pCV, pDV, pEV, pFV
1568      REAL(wp), DIMENSION(:,:), INTENT(inout) ::  pvel_res                       ! velocity residual @last inner iteration
1569      REAL(wp), DIMENSION(:,:), INTENT(inout) ::  pvel_diff                      ! velocity difference @last outer iteration
1570      !!
1571
1572      INTEGER           ::   idtime, istatus, ix_dim, iy_dim
1573      INTEGER           ::   ji, jj          ! dummy loop indices
1574      INTEGER           ::   it_inn_file, it_out_file
1575      REAL(wp)          ::   zu_res_mean, zv_res_mean, zvel_res_mean                  ! mean residuals of the linear system
1576      REAL(wp)          ::   zu_mad, zv_mad, zvel_mad                                 ! mean absolute deviation, sub-iterates
1577      REAL(wp)          ::   zu_mad_outer, zv_mad_outer, zvel_mad_outer               ! mean absolute deviation, outer-iterates
1578      REAL(wp)          ::   zvel_err_max, zmke, zu, zv                               ! local scalars
1579      REAL(wp)          ::   z1_pglob_area                                            ! inverse global ice area
1580
1581      REAL(wp), DIMENSION(jpi,jpj) ::   zu_res, zv_res, zvel2                         ! local arrays
1582      REAL(wp), DIMENSION(jpi,jpj) ::   zu_diff, zv_diff                              ! local arrays
1583                                                                             
1584      CHARACTER(len=20) ::   clname
1585      !!----------------------------------------------------------------------
1586
1587
1588      IF( lwp ) THEN
1589
1590         WRITE(numout,*)
1591         WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control'
1592         WRITE(numout,*) '~~~~~~~~~~~'
1593         WRITE(numout,*) ' kt          =  : ', kt
1594         WRITE(numout,*) ' kitout      =  : ', kitout
1595         WRITE(numout,*) ' kitinn      =  : ', kitinn
1596         WRITE(numout,*) ' kitinntot   =  : ', kitinntot
1597         WRITE(numout,*) ' kitoutmax (nn_nout_vp) =  ', kitoutmax
1598         WRITE(numout,*) ' kitinnmax (nn_ninn_vp) =  ', kitinnmax
1599         WRITE(numout,*) ' kitinntotmax (nn_nvp)  =  ', kitinntotmax
1600         WRITE(numout,*)
1601
1602      ENDIF
1603
1604      z1_pglob_area = 1._wp / pglob_area      ! inverse global ice area
1605
1606      ! create file
1607      IF( kt == nit000 .AND. kitinntot == 1 ) THEN
1608         !
1609         IF( lwm ) THEN
1610
1611            clname = 'ice_cvg.nc'
1612            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
1613            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
1614
1615            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
1616            istatus = NF90_DEF_DIM( ncvgid, 'x'     , jpi, ix_dim )
1617            istatus = NF90_DEF_DIM( ncvgid, 'y'     , jpj, iy_dim )
1618
1619            istatus = NF90_DEF_VAR( ncvgid, 'u_res'         , NF90_DOUBLE  , (/ idtime /), nvarid_ures )
1620            istatus = NF90_DEF_VAR( ncvgid, 'v_res'         , NF90_DOUBLE  , (/ idtime /), nvarid_vres )
1621            istatus = NF90_DEF_VAR( ncvgid, 'vel_res'       , NF90_DOUBLE  , (/ idtime /), nvarid_velres )
1622
1623            istatus = NF90_DEF_VAR( ncvgid, 'uerr_max_sub'  , NF90_DOUBLE  , (/ idtime /), nvarid_uerr_max )
1624            istatus = NF90_DEF_VAR( ncvgid, 'verr_max_sub'  , NF90_DOUBLE  , (/ idtime /), nvarid_verr_max )
1625            istatus = NF90_DEF_VAR( ncvgid, 'velerr_max_sub', NF90_DOUBLE  , (/ idtime /), nvarid_velerr_max )
1626
1627            istatus = NF90_DEF_VAR( ncvgid, 'umad_sub'      , NF90_DOUBLE  , (/ idtime /), nvarid_umad )
1628            istatus = NF90_DEF_VAR( ncvgid, 'vmad_sub'      , NF90_DOUBLE  , (/ idtime /), nvarid_vmad )
1629            istatus = NF90_DEF_VAR( ncvgid, 'velmad_sub'    , NF90_DOUBLE  , (/ idtime /), nvarid_velmad )
1630           
1631            istatus = NF90_DEF_VAR( ncvgid, 'umad_outer'    , NF90_DOUBLE  , (/ idtime /), nvarid_umad_outer   )
1632            istatus = NF90_DEF_VAR( ncvgid, 'vmad_outer'    , NF90_DOUBLE  , (/ idtime /), nvarid_vmad_outer   )
1633            istatus = NF90_DEF_VAR( ncvgid, 'velmad_outer'  , NF90_DOUBLE  , (/ idtime /), nvarid_velmad_outer )
1634
1635            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE  , (/ idtime /), nvarid_mke )
1636
1637            istatus = NF90_ENDDEF(ncvgid)
1638
1639         ENDIF
1640         !
1641      ENDIF
1642
1643      !------------------------------------------------------------
1644      !
1645      ! Max absolute velocity difference with previous sub-iterate
1646      ! ( zvel_err_max )
1647      !
1648      !------------------------------------------------------------
1649      !
1650      ! This comes from the criterion used to stop the iterative procedure
1651      zvel_err_max   = 0.5_wp * ( puerr_max + pverr_max ) ! average of U- and V- maximum error over the whole domain
1652
1653      !----------------------------------------------
1654      !
1655      ! Mean-absolute-deviation (sub-iterates)
1656      ! ( zu_mad, zv_mad, zvel_mad)
1657      !
1658      !----------------------------------------------
1659      !
1660      ! U
1661      zu_diff(:,:) = 0._wp
1662      DO jj = 2, jpj - 1
1663         DO ji = 2, jpi - 1
1664            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
1665         END DO
1666      END DO
1667
1668      ! V
1669      zv_diff(:,:)   = 0._wp
1670      DO jj = 2, jpj - 1
1671         DO ji = 2, jpi - 1
1672            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
1673         END DO
1674      END DO
1675
1676
1677      ! global sum & U-V average
1678      CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff,  'U',  1., zv_diff , 'V',  1. )
1679      zu_mad   = glob_sum( 'icedyn_rhg_vp : ', zu_diff )
1680      zv_mad   = glob_sum( 'icedyn_rhg_vp : ', zv_diff )
1681
1682      zvel_mad = 0.5_wp * ( zu_mad + zv_mad )
1683
1684      !-----------------------------------------------
1685      !
1686      ! Mean-absolute-deviation (outer-iterates)
1687      ! ( zu_mad_outer, zv_mad_outer, zvel_mad_outer)
1688      !
1689      !-----------------------------------------------
1690      !
1691      IF ( kitinn == kitinnmax ) THEN ! only work at the end of outer iterates
1692
1693         ! * U
1694         zu_diff(:,:) = 0._wp
1695         DO jj = 2, jpj - 1
1696            DO ji = 2, jpi - 1
1697               zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub_outer(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * &
1698                              &    z1_pglob_area
1699            END DO
1700         END DO
1701
1702         ! * V
1703         zv_diff(:,:)   = 0._wp
1704         DO jj = 2, jpj - 1
1705            DO ji = 2, jpi - 1
1706               zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb_outer(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * &
1707                              &    z1_pglob_area
1708            END DO
1709         END DO
1710
1711         ! Global ice-concentration, grid-cell-area weighted mean
1712         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff,  'U',  1., zv_diff , 'V',  1. ) ! abs behaves as a scalar no ?
1713
1714         zu_mad_outer   = glob_sum( 'icedyn_rhg_vp : ', zu_diff )
1715         zv_mad_outer   = glob_sum( 'icedyn_rhg_vp : ', zv_diff )
1716   
1717         ! Average of both U & V
1718         zvel_mad_outer = 0.5_wp * ( zu_mad_outer + zv_mad_outer )
1719                 
1720      ENDIF
1721
1722      ! --- Spatially-resolved absolute difference to send back to main routine
1723      ! (last iteration only, T-point)
1724
1725      IF ( kitinntot == kitinntotmax) THEN
1726
1727         zu_diff(:,:) = 0._wp
1728         zv_diff(:,:) = 0._wp
1729
1730         DO jj = 2, jpj - 1
1731            DO ji = 2, jpi - 1
1732
1733               zu_diff(ji,jj) = ( ABS ( ( pu(ji-1,jj) - pub_outer(ji-1,jj) ) ) * umask(ji-1,jj,1) &
1734    &                           + ABS ( ( pu(ji,jj  ) - pub_outer(ji,jj)   ) ) * umask(ji,jj,1) ) &
1735    &                           / ( umask(ji-1,jj,1) + umask(ji,jj,1) )
1736
1737               zv_diff(ji,jj) = ( ABS ( ( pv(ji,jj-1)   - pvb_outer(ji,jj-1) ) ) * vmask(ji,jj-1,1) &
1738    &                           + ABS ( ( pv(ji,jj  ) - pvb_outer(ji,jj)     ) ) * vmask(ji,jj,1)   &
1739    &                           / ( vmask(ji,jj-1,1) + vmask(ji,jj,1) ) )
1740     
1741            END DO
1742         END DO
1743
1744         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff,  'T',  1., zv_diff , 'T',  1. )
1745         pvel_diff(:,:) = 0.5_wp * ( zu_diff(:,:) + zv_diff(:,:) )
1746
1747      ELSE
1748
1749         pvel_diff(:,:) = 0._wp
1750
1751      ENDIF
1752
1753      !---------------------------------------
1754      !
1755      ! ---  Mean residual & kinetic energy
1756      !
1757      !---------------------------------------
1758
1759      IF ( kitinntot == 1 ) THEN
1760
1761         zu_res_mean   = 0._wp
1762         zv_res_mean   = 0._wp
1763         zvel_res_mean = 0._wp
1764         zmke          = 0._wp
1765
1766      ELSE
1767
1768         ! * Mean residual (N/m2)
1769         ! Here we take the residual of the linear system (N/m2),
1770         ! We define it as in mitgcm: global area-weighted mean of square-root residual
1771         ! Local residual r = Ax - B expresses to which extent the momentum balance is verified
1772         ! i.e., how close we are to a solution
1773         !!!! UNITS OF RESIDUAL ARE NOT m/s BUT N/m2 (CHECK)
1774         zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp
1775
1776         DO jj = 2, jpj - 1
1777            DO ji = 2, jpi - 1                                     
1778               zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               &
1779                  &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
1780               zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               &
1781                  &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
1782
1783!              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)
1784!              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)
1785   
1786               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
1787               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
1788   
1789            END DO
1790         END DO                 
1791
1792         ! Global ice-concentration, grid-cell-area weighted mean
1793         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_res,  'U', 1., zv_res , 'V', 1. )
1794   
1795         zu_res_mean   = glob_sum( 'ice_rhg_vp', zu_res(:,:) )
1796         zv_res_mean   = glob_sum( 'ice_rhg_vp', zv_res(:,:) )
1797         zvel_res_mean = 0.5_wp * ( zu_res_mean + zv_res_mean )
1798   
1799         ! --- Global mean kinetic energy per unit area (J/m2)
1800         zvel2(:,:) = 0._wp
1801         DO jj = 2, jpj - 1
1802            DO ji = 2, jpi - 1                   
1803               zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point
1804               zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) )
1805               zvel2(ji,jj)  = zu*zu + zv*zv              ! square of ice velocity at T-point 
1806            END DO
1807         END DO
1808         
1809         zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area
1810   
1811      ENDIF ! kitinntot
1812
1813      !--- Spatially-resolved residual at last iteration to send back to main routine (last iteration only)
1814      !--- Calculation @T-point
1815
1816      IF ( kitinntot == kitinntotmax) THEN
1817
1818         DO jj = 2, jpj - 1 ! @ residuals at u and v points
1819            DO ji = 2, jpi - 1
1820
1821               zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               &
1822                  &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
1823               zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               &
1824                  &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
1825
1826               zu_res(ji,jj)  = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) 
1827               zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) 
1828
1829            END DO
1830         END DO
1831
1832         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_res,  'U',  1., zv_res , 'V',  1. )
1833
1834         DO jj = 2, jpj - 1         ! average @T-point
1835            DO ji = 2, jpi - 1
1836               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) )
1837            END DO
1838         END DO
1839
1840         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. )
1841
1842      ELSE
1843
1844         pvel_res(:,:) = 0._wp
1845
1846      ENDIF
1847                 
1848      !                                                ! ==================== !
1849
1850      it_inn_file =  ( kt - 1 ) * kitinntotmax + kitinntot ! time step in the file
1851      it_out_file =  ( kt - 1 ) * kitoutmax  + kitout
1852
1853      ! write variables
1854      IF( lwm ) THEN
1855
1856         istatus = NF90_PUT_VAR( ncvgid, nvarid_ures  , (/zu_res_mean/), (/it_inn_file/), (/1/) )        ! Residuals of the linear system, area weighted mean
1857         istatus = NF90_PUT_VAR( ncvgid, nvarid_vres  , (/zv_res_mean/), (/it_inn_file/), (/1/) )        !
1858         istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvel_res_mean/), (/it_inn_file/), (/1/) )      !
1859
1860         istatus = NF90_PUT_VAR( ncvgid, nvarid_uerr_max  , (/puerr_max/), (/it_inn_file/), (/1/) )      ! Max velocit_inn_filey error, sub-it_inn_fileerates
1861         istatus = NF90_PUT_VAR( ncvgid, nvarid_verr_max  , (/pverr_max/), (/it_inn_file/), (/1/) )      !
1862         istatus = NF90_PUT_VAR( ncvgid, nvarid_velerr_max, (/zvel_err_max/), (/it_inn_file/), (/1/) )   !
1863
1864         istatus = NF90_PUT_VAR( ncvgid, nvarid_umad    , (/zu_mad/)  , (/it_inn_file/), (/1/) )         ! velocit_inn_filey MAD, area/sic-weighted, sub-it_inn_fileerates
1865         istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad    , (/zv_mad/)  , (/it_inn_file/), (/1/) )         !
1866         istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad  , (/zvel_mad/), (/it_inn_file/), (/1/) )         !
1867
1868         istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/kitinntot/), (/1/) )                    ! mean kinetic energy
1869
1870         IF ( kitinn == kitinnmax ) THEN ! only print outer mad at the end of inner loop
1871
1872            istatus = NF90_PUT_VAR( ncvgid, nvarid_umad_outer    , (/zu_mad_outer/)  , (/it_out_file/), (/1/) )   ! velocity MAD, area/sic-weighted, outer-iterates
1873            istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad_outer    , (/zv_mad_outer/)  , (/it_out_file/), (/1/) )   !
1874            istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad_outer  , (/zvel_mad_outer/), (/it_out_file/), (/1/) )   !
1875
1876         ENDIF
1877
1878         ! spatially-dependent things (to be coded)
1879         ! velocity difference ?
1880         ! residual ?
1881!        !
1882!        IF ( kiter == kitermax ) THEN
1883!           WRITE(numout,*) ' Should plot the spatially dependent residual '
1884!           istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) )          ! U-residual, spatially dependent
1885!           istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) )          ! V-residual, spatially dependent
1886!        ENDIF
1887
1888         ! close file
1889         IF( kt == nitend )   istatus = NF90_CLOSE( ncvgid )
1890      ENDIF
1891     
1892   END SUBROUTINE rhg_cvg_vp
1893   
1894
1895   
1896#else
1897   !!----------------------------------------------------------------------
1898   !!   Default option         Empty module           NO SI3 sea-ice model
1899   !!----------------------------------------------------------------------
1900#endif
1901
1902   !!==============================================================================
1903END MODULE icedyn_rhg_vp
1904
Note: See TracBrowser for help on using the repository browser.