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

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

Further fixes to 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. )
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. )
908 
909                     !-----------------------------
910                     ! Backward sweep & relaxation
911                     !-----------------------------
912
913                     DO jj = jj_min, jpj - 1, nn_zebra_vp
914                   
915                        ! --- Backward sweep
916                        ! last row
917                        zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) )
918                        u_ice(jpi-1,jj)    = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 
919                                           &            * umask(jpi-1,jj,1)
920                        DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop
921                           zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) )
922                           u_ice(ji,jj)    = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1)   & 
923                                           &                  / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 ) 
924                        END DO
925
926                        !--- Relaxation
927                        ! and velocity masking for little-ice and no-ice cases
928                        DO ji = 2, jpi - 1   
929                       
930                           u_ice(ji,jj) = zu_b(ji,jj) + rn_relaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) ! relaxation
931                           
932                           u_ice(ji,jj) =   zmsk00x(ji,jj)                                        &   ! masking
933                              &         * (           zmsk01x(ji,jj)   * u_ice(ji,jj)             &
934                              &           + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp     ) * umask(ji,jj,1)
935                           
936                        END DO
937
938                     END DO ! jj
939                     
940                  END DO ! zebra loop
941
942               ENDIF !   ll_u_iterate
943
944               !                           ! ---------------------------- !
945               IF ( ll_v_iterate ) THEN    ! --- Solve for V-velocity --- !
946               !                           ! ---------------------------- !
947                                         
948                  ! MV OPT: what follows could be subroutinized...                 
949                  ! Thomas Algorithm  for tridiagonal solver
950                  ! A*v(i,j-1)+B*v(i,j)+C*v(i,j+1) = F
951                  ! It is intentional to have a ji then jj loop for V-velocity
952                  !!! ZH97 explain it is critical for convergence speed
953
954                  zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; zCV_prime(:,:) = 0._wp
955
956                  DO jn = 1, nn_zebra_vp ! "zebra" loop
957     
958                     IF ( jn == 1 ) THEN   ;   ji_min = 2 
959                     ELSE                  ;   ji_min = 3
960                     ENDIF
961
962                     IF ( lwp ) WRITE(numout,*) ' Into the V-zebra loop at step jn = ', jn, ', with ji_min = ', ji_min
963         
964                     DO ji = ji_min, jpi - 1, nn_zebra_vp 
965                     
966                        !------------------------
967                        ! Independent term (zFV)
968                        !------------------------
969                        DO jj = 2, jpj - 1
970
971                           ! boundary condition substitution (check it is correctly applied !!!)
972                           ! see Zhang and Hibler, 1997, Appendix B
973                           zAA3 = 0._wp
974                           IF ( jj == 2 )       zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1)
975                           IF ( jj == jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1)
976     
977                           ! right hand side
978                           zFV(ji,jj) = ( zrhsv(ji,jj) &                                   ! right-hand side terms
979                               &        + zAA3                                           & ! boundary condition translation
980                               &        + zDV(ji,jj) * v_ice(ji-1,jj)                    & ! internal force, j-1
981                               &        + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1)   ! internal force, j+1
982
983                        END DO
984                       
985                     END DO
986
987                     CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV, 'V',  1.)
988                     
989                     !---------------
990                     ! Forward sweep
991                     !---------------
992                     DO ji = ji_min, jpi - 1, nn_zebra_vp 
993                     
994                        ! MV BUG OCT 6.. Both zBV_prime and zFV_prime were undefined for jj = 2
995                        zBV_prime(ji,2)     = zBV(ji,2)
996                        zFV_prime(ji,2)     = zFV(ji,2)
997
998                        DO jj = 3, jpj - 1 
999
1000                           zfac             = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) )
1001                           zw               = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 )
1002                           zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1)
1003                           zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 
1004
1005                        END DO
1006
1007                     END DO
1008
1009                     CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV_prime, 'V',  1., zBV_prime, 'V', 1. )
1010                     
1011                     !-----------------------------
1012                     ! Backward sweep & relaxation
1013                     !-----------------------------
1014                     DO ji = ji_min, jpi - 1, nn_zebra_vp 
1015                   
1016                        ! --- Backward sweep
1017                        ! last row
1018                        zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) )
1019                        v_ice(ji,jpj-1)  = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) & 
1020                                         &         * vmask(ji,jpj-1,1)  ! last row
1021
1022                        ! other rows
1023                        DO jj = jpj-2, 2, -1 ! original back loop
1024                           zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) )
1025                           v_ice(ji,jj)   = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) &
1026                                          &  / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 )       
1027                        END DO           
1028                                                   
1029                        ! --- Relaxation  & masking (should it be now or later)
1030                        DO jj = 2, jpj - 1
1031                       
1032                            v_ice(ji,jj) = zv_b(ji,jj) + rn_relaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) )    ! relaxation
1033                           
1034                            v_ice(ji,jj) =   zmsk00y(ji,jj)                                        &      ! masking
1035                              &         * (           zmsk01y(ji,jj)   * v_ice(ji,jj)              &
1036                              &           + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) * vmask(ji,jj,1)
1037
1038                        END DO ! jj
1039                       
1040                     END DO ! ji
1041                     
1042                  END DO ! zebra loop
1043                                   
1044               ENDIF !   ll_v_iterate
1045
1046               CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
1047                             
1048               !--------------------------------------------------------------------------------------
1049               ! -- Check convergence based on maximum velocity difference, continue or stop the loop
1050               !--------------------------------------------------------------------------------------
1051
1052               !------
1053               ! on U
1054               !------
1055               ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check
1056               ! then we must optimize the use of the mpp_max, which is prohibitive                           
1057               zuerr_max  = 0._wp
1058                               
1059               IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN
1060
1061                  ! - Maximum U-velocity difference               
1062                  zuerr(:,:) = 0._wp
1063                  DO jj = 2, jpj - 1
1064                     DO ji = 2, jpi - 1
1065                        zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 
1066                     END DO
1067                  END DO
1068                  zuerr_max = MAXVAL( zuerr )
1069                  CALL mpp_max( 'icedyn_rhg_evp', zuerr_max )   ! max over the global domain - damned!
1070
1071                  ! - Stop if max :error is too large ("safeguard against bad forcing" of original Zhang routine)
1072                  IF ( i_inn > 1 .AND. zuerr_max > rn_uerr_max_vp ) THEN
1073                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
1074                      ll_u_iterate = .FALSE.
1075                  ENDIF
1076                 
1077                  ! - Stop if error small enough
1078                  IF ( zuerr_max < rn_uerr_min_vp ) THEN                                       
1079                      IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! "
1080                      ll_u_iterate = .FALSE.
1081                  ENDIF
1082                                               
1083               ENDIF ! ll_u_iterate
1084
1085               !------
1086               ! on V
1087               !------
1088               zverr_max = 0._wp
1089               
1090               IF ( ll_v_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN
1091               
1092                  ! - Maximum V-velocity difference
1093                  zverr(:,:)   = 0._wp   
1094                  DO jj = 2, jpj - 1
1095                     DO ji = 2, jpi - 1
1096                        zverr(ji,jj) = ABS ( ( v_ice(ji,jj) - zv_b(ji,jj) ) ) * vmask(ji,jj,1)
1097                     END DO
1098                  END DO
1099                  zverr_max = MAXVAL( zverr )
1100                  CALL mpp_max( 'icedyn_rhg_evp', zverr_max )   ! max over the global domain - damned!
1101
1102                  ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)
1103                  IF ( i_inn > 1 .AND. zverr_max > rn_uerr_max_vp ) THEN
1104                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zverr_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
1105                      ll_v_iterate = .FALSE.
1106                  ENDIF
1107                 
1108                  ! - Stop if error small enough
1109                  IF ( zverr_max < rn_uerr_min_vp ) THEN                                       
1110                      IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! "
1111                      ll_v_iterate = .FALSE.
1112                  ENDIF
1113                 
1114               ENDIF ! ll_v_iterate
1115
1116            ENDIF ! ---    end ll_u_iterate or ll_v_iterate
1117               
1118            !---------------------------------------------------------------------------------------
1119            !
1120            ! --- Calculate extra convergence diagnostics and save them
1121            !
1122            !---------------------------------------------------------------------------------------
1123            IF( ln_rhg_chkcvg .AND. MOD ( i_inn - 1, nn_cvgchk_vp ) == 0 ) THEN
1124
1125               CALL rhg_cvg_vp( kt, i_out, i_inn, i_inn_tot, nn_nout_vp, nn_ninn_vp, nn_nvp,        &
1126                      &         u_ice, v_ice, zu_b, zv_b, zu_b_outer, zv_b_outer,                   &
1127                      &         zmt, za_iU, za_iV, zuerr_max, zverr_max, zglob_area,                &
1128                      &         zrhsu, zAU, zBU, zCU, zDU, zEU, zFU,                                &
1129                      &         zrhsv, zAV, zBV, zCV, zDV, zEV, zFV,                                &
1130                                zvel_res, zvel_diff )
1131
1132            ENDIF
1133
1134            IF ( lwp ) WRITE(numout,*) ' Done convergence tests '
1135
1136         END DO ! i_inn, end of inner loop
1137
1138      END DO ! End of outer loop (i_out) =============================================================================================
1139
1140      IF ( lwp ) WRITE(numout,*) ' We are out of outer loop '
1141
1142      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU  , 'U',  1., zFV  , 'V',  1. )
1143      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU_prime  , 'U',  1., zBV_prime  , 'V',  1. )
1144      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU_prime  , 'U',  1., zFV_prime  , 'V',  1. )
1145      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCU_prime  , 'U',  1., zCV_prime  , 'V',  1. )
1146
1147      CALL iom_put( 'zFU'           , zFU            ) ! MV DEBUG
1148      CALL iom_put( 'zBU_prime'     , zBU_prime      ) ! MV DEBUG
1149      CALL iom_put( 'zCU_prime'     , zCU_prime      ) ! MV DEBUG
1150      CALL iom_put( 'zFU_prime'     , zFU_prime      ) ! MV DEBUG
1151
1152      CALL iom_put( 'zFV'           , zFV            ) ! MV DEBUG
1153      CALL iom_put( 'zBV_prime'     , zBV_prime      ) ! MV DEBUG
1154      CALL iom_put( 'zCV_prime'     , zCV_prime      ) ! MV DEBUG
1155      CALL iom_put( 'zFV_prime'     , zFV_prime      ) ! MV DEBUG
1156
1157      CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
1158
1159
1160      IF( ln_rhg_chkcvg ) THEN
1161         
1162         IF( iom_use('velo_res') )   CALL iom_put( 'velo_res', zvel_res  )   ! linear system residual  @last inner&outer iteration
1163         IF( iom_use('velo_ero') )   CALL iom_put( 'velo_ero', zvel_diff )   ! abs velocity difference @last outer iteration
1164         IF( iom_use('uice_eri') )   CALL iom_put( 'uice_eri', zuerr     )   ! abs velocity difference @last inner iteration
1165         IF( iom_use('vice_eri') )   CALL iom_put( 'vice_eri', zverr     )   ! abs velocity difference @last inner iteration
1166
1167         DEALLOCATE( zvel_res , zvel_diff )
1168       
1169      ENDIF ! ln_rhg_chkcvg
1170
1171      ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead
1172      ! info: probably wrong because when doing du/dt = div ( sigma ) we have spurious waves instead of u = v = 0
1173
1174!     DO jj = 2, jpj - 1
1175!        DO ji = 2, jpi - 1   
1176
1177!            u_ice(ji,jj) =   zmsk00x(ji,jj)                               &
1178!     &         * (           zmsk01x(ji,jj)   * u_oce(ji,jj) * 0.01_wp    &
1179!                 + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    )
1180
1181!            v_ice(ji,jj) =   zmsk00y(ji,jj)                               &
1182!     &         * (           zmsk01y(ji,jj)   * v_oce(ji,jj) * 0.01_wp    &
1183!                 + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    )
1184
1185!        END DO
1186!     END DO
1187
1188!     CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
1189
1190!     IF ( lwp ) WRITE(numout,*) ' Velocity replaced '
1191
1192      ! END DEBUG
1193
1194      !------------------------------------------------------------------------------!
1195      !
1196      ! --- Recompute delta, shear and div (inputs for mechanical redistribution)
1197      !
1198      !------------------------------------------------------------------------------!
1199      !
1200      ! MV OPT: subroutinize ?
1201
1202      DO jj = 1, jpj - 1
1203         DO ji = 1, jpi - 1
1204
1205            ! shear at F points
1206            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)   &
1207               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
1208               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
1209
1210         END DO
1211      END DO           
1212     
1213      DO jj = 2, jpj - 1
1214         DO ji = 2, jpi - 1 !
1215           
1216            ! tension**2 at T points
1217            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)   &
1218               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
1219               &   ) * r1_e1e2t(ji,jj)
1220            zdt2 = zdt * zdt
1221           
1222            zten_i(ji,jj) = zdt
1223           
1224            ! shear**2 at T points (doc eq. A16)
1225            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
1226               &   + 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)  &
1227               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
1228           
1229            ! shear at T points
1230            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
1231
1232            ! divergence at T points
1233            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
1234               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
1235               &             ) * r1_e1e2t(ji,jj)
1236           
1237            ! delta at T points
1238            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
1239            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
1240           
1241            !pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
1242            pdelta_i(ji,jj) = zdelta + rn_creepl
1243
1244         END DO
1245      END DO
1246
1247      IF ( lwp ) WRITE(numout,*) ' Deformation recalculated '
1248     
1249      CALL lbc_lnk_multi( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
1250     
1251      !------------------------------------------------------------------------------!
1252      !
1253      ! --- Diagnostics
1254      !
1255      !------------------------------------------------------------------------------!
1256      !
1257      ! MV OPT: subroutinize ?
1258      !
1259      !----------------------------------
1260      ! --- Recompute stresses if needed
1261      !----------------------------------
1262      !
1263      ! ---- Sea ice stresses at T-points
1264      IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1265     
1266         DO jj = 2, jpj - 1
1267            DO ji = 2, jpi - 1
1268               zp_deltastar_t(ji,jj)   =   strength(ji,jj) / pdelta_i(ji,jj) 
1269               zfac                    =   zp_deltastar_t(ji,jj) 
1270               zs1(ji,jj)              =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
1271               zs2(ji,jj)              =   zfac * z1_ecc2 * zten_i(ji,jj)
1272               zs12(ji,jj)             =   zfac * z1_ecc2 * pshear_i(ji,jj)
1273            END DO
1274         END DO
1275
1276         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. )
1277     
1278      ENDIF
1279     
1280      ! ---- s12 at F-points     
1281      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1282
1283         DO jj = 1, jpj - 1
1284            DO ji = 1, jpi - 1
1285           
1286               ! P/delta* at F points
1287               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) )
1288               
1289               ! s12 at F-points
1290               zs12f(ji,jj) = zp_deltastar_f * z1_ecc2 * zds(ji,jj)
1291               
1292            END DO
1293         END DO
1294
1295         CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. )
1296         
1297      ENDIF
1298
1299      IF ( lwp ) WRITE(numout,*) ' zs12f recalculated '
1300
1301      !
1302      !-----------------------
1303      ! --- Store diagnostics
1304      !-----------------------
1305      !
1306      ! --- Ice-ocean, ice-atm. & ice-ocean bottom (landfast) stresses --- !
1307      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
1308         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
1309
1310         ALLOCATE( ztaux_oi(jpi,jpj) , ztauy_oi(jpi,jpj) )
1311
1312         !--- Recalculate oceanic stress at last inner iteration
1313         DO jj = 2, jpj - 1
1314            DO ji = 2, jpi - 1
1315
1316                !--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
1317                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)
1318                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)
1319               
1320                !--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
1321                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) )  &
1322                  &                                                + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
1323                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) )  &
1324                  &                                                + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
1325                 
1326                !--- Ocean-ice stress
1327                ztaux_oi(ji,jj) = zCwU(ji,jj) * ( u_oce(ji,jj) - u_ice(ji,jj) )
1328                ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) )
1329               
1330            END DO
1331         END DO
1332         
1333         !
1334         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, &
1335!            &                                 ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
1336         !
1337         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
1338         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
1339         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
1340         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
1341!        CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
1342!        CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
1343
1344         DEALLOCATE( ztaux_oi , ztauy_oi )
1345
1346      ENDIF
1347       
1348      ! --- Divergence, shear and strength --- !
1349      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
1350      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
1351      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
1352      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
1353
1354      IF ( lwp ) WRITE(numout,*) 'Some terms recalculated '
1355
1356      ! --- Stress tensor invariants (SIMIP diags) --- !
1357      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
1358         !
1359         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags.
1360         ! Definitions following Coon (1974) and Feltham (2008)
1361         !
1362         ! sigma1, sigma2, sigma12 are useful (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
1363         ! however these are NOT stress tensor components, neither stress invariants, nor stress principal components
1364         !
1365         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
1366         !         
1367         DO jj = 2, jpj - 1
1368            DO ji = 2, jpi - 1
1369               ! Stress invariants
1370               zsig_I(ji,jj)    =   zs1(ji,jj) * 0.5_wp                                        ! 1st invariant, aka average normal stress aka negative pressure
1371               zsig_II(ji,jj)   =   SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj) )  ! 2nd invariant, aka maximum shear stress               
1372            END DO
1373         END DO
1374
1375         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
1376         
1377         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,   zsig_I(:,:)  * zmsk00(:,:) ) ! Normal stress
1378         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' ,   zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
1379         
1380         DEALLOCATE ( zsig_I, zsig_II )
1381         
1382      ENDIF
1383
1384      IF ( lwp ) WRITE(numout,*) 'SIMIP work done'
1385
1386      ! --- Normalized stress tensor principal components --- !
1387      ! These are used to plot the normalized yield curve (Lemieux & Dupont, GMD 2020)
1388      ! To plot the yield curve and evaluate physical convergence, they have two recommendations
1389      ! Recommendation 1 : Use ice strength, not replacement pressure
1390      ! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765)
1391      ! R2 means we need to recompute stresses
1392
1393      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
1394         !
1395         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
1396         !         
1397         DO jj = 2, jpj - 1
1398            DO ji = 2, jpi - 1
1399               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
1400               !                        and **deformations** at current iterates
1401               !                        following Lemieux & Dupont (2020)
1402               zfac             =   zp_deltastar_t(ji,jj)
1403               zsig1            =   zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) )
1404               zsig1            = 0._wp !!! FUCKING DEBUG TEST !!!
1405               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
1406               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
1407               
1408               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
1409               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                              ! 1st invariant
1410               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd invariant
1411
1412               ! Normalized  principal stresses (used to display the ellipse)
1413               z1_strength      =   1._wp / MAX ( 1._wp , strength(ji,jj) )
1414               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
1415               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
1416            END DO
1417         END DO
1418         IF ( lwp ) WRITE(numout,*) 'Some shitty stress work done'
1419         !
1420         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
1421         !     
1422         IF ( lwp ) WRITE(numout,*) ' Beauaaaarflblbllll '
1423         !
1424         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
1425         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
1426
1427         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II )
1428
1429         IF ( lwp ) WRITE(numout,*) ' So what ??? '
1430         
1431      ENDIF
1432
1433      ! --- SIMIP, terms of tendency for momentum equation  --- !
1434      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
1435         & iom_use('corstrx') .OR. iom_use('corstry') ) THEN
1436
1437         ! --- Recalculate Coriolis stress at last inner iteration
1438         DO jj = 2, jpj - 1
1439            DO ji = 2, jpi - 1
1440                ! --- U-component
1441                zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  &
1442                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
1443                           &             + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
1444                zCorV(ji,jj)         = - 0.25_wp * r1_e2v(ji,jj) *  &
1445                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
1446                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
1447            END DO
1448         END DO
1449         !
1450         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., &
1451            &                                 zCorU, 'U', -1., zCorV, 'V', -1. )
1452         !
1453         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
1454         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
1455         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
1456         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
1457
1458      ENDIF
1459     
1460      IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1461
1462         ! Recalculate internal forces (divergence of stress tensor) at last inner iteration
1463         DO jj = 2, jpj - 1
1464            DO ji = 2, jpi - 1
1465
1466               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
1467                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
1468                  &                    ) * r1_e2u(ji,jj)                                                                      &
1469                  &                  + ( zs12f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
1470                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
1471                  &                  ) * r1_e1e2u(ji,jj)
1472
1473               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
1474                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
1475                  &                    ) * r1_e1v(ji,jj)                                                                      &
1476                  &                  + ( zs12f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
1477                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
1478                  &                  ) * r1_e1e2v(ji,jj)
1479
1480            END DO
1481         END DO
1482           
1483         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. )
1484         
1485         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
1486         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
1487         
1488      ENDIF
1489
1490      ! --- Ice & snow mass and ice area transports
1491      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
1492         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
1493         !
1494         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
1495            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
1496         !
1497         DO jj = 2, jpj - 1
1498            DO ji = 2, jpi - 1
1499               ! 2D ice mass, snow mass, area transport arrays (X, Y)
1500               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
1501               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
1502
1503               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
1504               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
1505
1506               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
1507               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
1508
1509               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
1510               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
1511
1512            END DO
1513         END DO
1514
1515         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
1516            &                                 zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
1517            &                                 zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
1518
1519         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
1520         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport
1521         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
1522         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
1523         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
1524         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
1525
1526         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
1527            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
1528
1529      ENDIF
1530
1531      DEALLOCATE( zmsk00, zmsk15 )
1532
1533   END SUBROUTINE ice_dyn_rhg_vp
1534   
1535   
1536   SUBROUTINE rhg_cvg_vp( kt, kitout, kitinn, kitinntot, kitoutmax, kitinnmax, kitinntotmax , &
1537                  &       pu, pv, pub, pvb, pub_outer, pvb_outer                     , &
1538                  &       pmt, pat_iu, pat_iv, puerr_max, pverr_max, pglob_area      , &
1539                  &       prhsu, pAU, pBU, pCU, pDU, pEU, pFU                        , &
1540                  &       prhsv, pAV, pBV, pCV, pDV, pEV, pFV                        , &   
1541                  &       pvel_res, pvel_diff                                            )
1542      !!
1543      !!----------------------------------------------------------------------
1544      !!                    ***  ROUTINE rhg_cvg_vp  ***
1545      !!                     
1546      !! ** Purpose :   check convergence of VP ice rheology
1547      !!
1548      !! ** Method  :   create a file ice_cvg.nc containing a few convergence diagnostics
1549      !!                This routine is called every sub-iteration, so it is cpu expensive
1550      !!
1551      !!                Calculates / stores
1552      !!                   - maximum absolute U-V difference (uice_cvg, u_dif, v_dif, m/s)
1553      !!                   - 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)
1554      !!                   - mean kinetic energy (mke_ice, J/m2)
1555      !!
1556      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)   
1557      !!
1558      !!----------------------------------------------------------------------
1559      !!
1560      INTEGER ,                 INTENT(in) ::   kt, kitout, kitinn, kitinntot    ! ocean model iterate, outer, inner and total n-iterations
1561      INTEGER ,                 INTENT(in) ::   kitoutmax, kitinnmax             ! max number of outer & inner iterations
1562      INTEGER ,                 INTENT(in) ::   kitinntotmax                     ! max number of total sub-iterations
1563      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb                 ! now & sub-iter-before velocities
1564      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pub_outer, pvb_outer             ! velocities @before outer iterations
1565      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pmt, pat_iu, pat_iv              ! mass at T-point, ice concentration at U&V
1566      REAL(wp),                 INTENT(in) ::   puerr_max, pverr_max             ! absolute mean velocity difference
1567      REAL(wp),                 INTENT(in) ::   pglob_area                       ! global ice area
1568      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsu, pAU, pBU, pCU, pDU, pEU, pFU ! linear system coefficients
1569      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsv, pAV, pBV, pCV, pDV, pEV, pFV
1570      REAL(wp), DIMENSION(:,:), INTENT(inout) ::  pvel_res                       ! velocity residual @last inner iteration
1571      REAL(wp), DIMENSION(:,:), INTENT(inout) ::  pvel_diff                      ! velocity difference @last outer iteration
1572      !!
1573
1574      INTEGER           ::   idtime, istatus, ix_dim, iy_dim
1575      INTEGER           ::   ji, jj          ! dummy loop indices
1576      INTEGER           ::   it_inn_file, it_out_file
1577      REAL(wp)          ::   zu_res_mean, zv_res_mean, zvel_res_mean                  ! mean residuals of the linear system
1578      REAL(wp)          ::   zu_mad, zv_mad, zvel_mad                                 ! mean absolute deviation, sub-iterates
1579      REAL(wp)          ::   zu_mad_outer, zv_mad_outer, zvel_mad_outer               ! mean absolute deviation, outer-iterates
1580      REAL(wp)          ::   zvel_err_max, zmke, zu, zv                               ! local scalars
1581      REAL(wp)          ::   z1_pglob_area                                            ! inverse global ice area
1582
1583      REAL(wp), DIMENSION(jpi,jpj) ::   zu_res, zv_res, zvel2                         ! local arrays
1584      REAL(wp), DIMENSION(jpi,jpj) ::   zu_diff, zv_diff                              ! local arrays
1585                                                                             
1586      CHARACTER(len=20) ::   clname
1587      !!----------------------------------------------------------------------
1588
1589
1590      IF( lwp ) THEN
1591
1592         WRITE(numout,*)
1593         WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control'
1594         WRITE(numout,*) '~~~~~~~~~~~'
1595         WRITE(numout,*) ' kt          =  : ', kt
1596         WRITE(numout,*) ' kitout      =  : ', kitout
1597         WRITE(numout,*) ' kitinn      =  : ', kitinn
1598         WRITE(numout,*) ' kitinntot   =  : ', kitinntot
1599         WRITE(numout,*) ' kitoutmax (nn_nout_vp) =  ', kitoutmax
1600         WRITE(numout,*) ' kitinnmax (nn_ninn_vp) =  ', kitinnmax
1601         WRITE(numout,*) ' kitinntotmax (nn_nvp)  =  ', kitinntotmax
1602         WRITE(numout,*)
1603
1604      ENDIF
1605
1606      z1_pglob_area = 1._wp / pglob_area      ! inverse global ice area
1607
1608      ! create file
1609      IF( kt == nit000 .AND. kitinntot == 1 ) THEN
1610         !
1611         IF( lwm ) THEN
1612
1613            clname = 'ice_cvg.nc'
1614            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
1615            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
1616
1617            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
1618            istatus = NF90_DEF_DIM( ncvgid, 'x'     , jpi, ix_dim )
1619            istatus = NF90_DEF_DIM( ncvgid, 'y'     , jpj, iy_dim )
1620
1621            istatus = NF90_DEF_VAR( ncvgid, 'u_res'         , NF90_DOUBLE  , (/ idtime /), nvarid_ures )
1622            istatus = NF90_DEF_VAR( ncvgid, 'v_res'         , NF90_DOUBLE  , (/ idtime /), nvarid_vres )
1623            istatus = NF90_DEF_VAR( ncvgid, 'vel_res'       , NF90_DOUBLE  , (/ idtime /), nvarid_velres )
1624
1625            istatus = NF90_DEF_VAR( ncvgid, 'uerr_max_sub'  , NF90_DOUBLE  , (/ idtime /), nvarid_uerr_max )
1626            istatus = NF90_DEF_VAR( ncvgid, 'verr_max_sub'  , NF90_DOUBLE  , (/ idtime /), nvarid_verr_max )
1627            istatus = NF90_DEF_VAR( ncvgid, 'velerr_max_sub', NF90_DOUBLE  , (/ idtime /), nvarid_velerr_max )
1628
1629            istatus = NF90_DEF_VAR( ncvgid, 'umad_sub'      , NF90_DOUBLE  , (/ idtime /), nvarid_umad )
1630            istatus = NF90_DEF_VAR( ncvgid, 'vmad_sub'      , NF90_DOUBLE  , (/ idtime /), nvarid_vmad )
1631            istatus = NF90_DEF_VAR( ncvgid, 'velmad_sub'    , NF90_DOUBLE  , (/ idtime /), nvarid_velmad )
1632           
1633            istatus = NF90_DEF_VAR( ncvgid, 'umad_outer'    , NF90_DOUBLE  , (/ idtime /), nvarid_umad_outer   )
1634            istatus = NF90_DEF_VAR( ncvgid, 'vmad_outer'    , NF90_DOUBLE  , (/ idtime /), nvarid_vmad_outer   )
1635            istatus = NF90_DEF_VAR( ncvgid, 'velmad_outer'  , NF90_DOUBLE  , (/ idtime /), nvarid_velmad_outer )
1636
1637            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE  , (/ idtime /), nvarid_mke )
1638
1639            istatus = NF90_ENDDEF(ncvgid)
1640
1641         ENDIF
1642         !
1643      ENDIF
1644
1645      !------------------------------------------------------------
1646      !
1647      ! Max absolute velocity difference with previous sub-iterate
1648      ! ( zvel_err_max )
1649      !
1650      !------------------------------------------------------------
1651      !
1652      ! This comes from the criterion used to stop the iterative procedure
1653      zvel_err_max   = 0.5_wp * ( puerr_max + pverr_max ) ! average of U- and V- maximum error over the whole domain
1654
1655      !----------------------------------------------
1656      !
1657      ! Mean-absolute-deviation (sub-iterates)
1658      ! ( zu_mad, zv_mad, zvel_mad)
1659      !
1660      !----------------------------------------------
1661      !
1662      ! U
1663      zu_diff(:,:) = 0._wp
1664      DO jj = 2, jpj - 1
1665         DO ji = 2, jpi - 1
1666            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
1667         END DO
1668      END DO
1669
1670      ! V
1671      zv_diff(:,:)   = 0._wp
1672      DO jj = 2, jpj - 1
1673         DO ji = 2, jpi - 1
1674            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
1675         END DO
1676      END DO
1677
1678
1679      ! global sum & U-V average
1680      CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff,  'U',  1., zv_diff , 'V',  1. )
1681      zu_mad   = glob_sum( 'icedyn_rhg_vp : ', zu_diff )
1682      zv_mad   = glob_sum( 'icedyn_rhg_vp : ', zv_diff )
1683
1684      zvel_mad = 0.5_wp * ( zu_mad + zv_mad )
1685
1686      !-----------------------------------------------
1687      !
1688      ! Mean-absolute-deviation (outer-iterates)
1689      ! ( zu_mad_outer, zv_mad_outer, zvel_mad_outer)
1690      !
1691      !-----------------------------------------------
1692      !
1693      IF ( kitinn == kitinnmax ) THEN ! only work at the end of outer iterates
1694
1695         ! * U
1696         zu_diff(:,:) = 0._wp
1697         DO jj = 2, jpj - 1
1698            DO ji = 2, jpi - 1
1699               zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub_outer(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * &
1700                              &    z1_pglob_area
1701            END DO
1702         END DO
1703
1704         ! * V
1705         zv_diff(:,:)   = 0._wp
1706         DO jj = 2, jpj - 1
1707            DO ji = 2, jpi - 1
1708               zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb_outer(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * &
1709                              &    z1_pglob_area
1710            END DO
1711         END DO
1712
1713         ! Global ice-concentration, grid-cell-area weighted mean
1714         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff,  'U',  1., zv_diff , 'V',  1. ) ! abs behaves as a scalar no ?
1715
1716         zu_mad_outer   = glob_sum( 'icedyn_rhg_vp : ', zu_diff )
1717         zv_mad_outer   = glob_sum( 'icedyn_rhg_vp : ', zv_diff )
1718   
1719         ! Average of both U & V
1720         zvel_mad_outer = 0.5_wp * ( zu_mad_outer + zv_mad_outer )
1721                 
1722      ENDIF
1723
1724      ! --- Spatially-resolved absolute difference to send back to main routine
1725      ! (last iteration only, T-point)
1726
1727      IF ( kitinntot == kitinntotmax) THEN
1728
1729         zu_diff(:,:) = 0._wp
1730         zv_diff(:,:) = 0._wp
1731
1732         DO jj = 2, jpj - 1
1733            DO ji = 2, jpi - 1
1734
1735               zu_diff(ji,jj) = ( ABS ( ( pu(ji-1,jj) - pub_outer(ji-1,jj) ) ) * umask(ji-1,jj,1) &
1736    &                           + ABS ( ( pu(ji,jj  ) - pub_outer(ji,jj)   ) ) * umask(ji,jj,1) ) &
1737    &                           / ( umask(ji-1,jj,1) + umask(ji,jj,1) )
1738
1739               zv_diff(ji,jj) = ( ABS ( ( pv(ji,jj-1)   - pvb_outer(ji,jj-1) ) ) * vmask(ji,jj-1,1) &
1740    &                           + ABS ( ( pv(ji,jj  ) - pvb_outer(ji,jj)     ) ) * vmask(ji,jj,1)   &
1741    &                           / ( vmask(ji,jj-1,1) + vmask(ji,jj,1) ) )
1742     
1743            END DO
1744         END DO
1745
1746         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff,  'T',  1., zv_diff , 'T',  1. )
1747         pvel_diff(:,:) = 0.5_wp * ( zu_diff(:,:) + zv_diff(:,:) )
1748
1749      ELSE
1750
1751         pvel_diff(:,:) = 0._wp
1752
1753      ENDIF
1754
1755      !---------------------------------------
1756      !
1757      ! ---  Mean residual & kinetic energy
1758      !
1759      !---------------------------------------
1760
1761      IF ( kitinntot == 1 ) THEN
1762
1763         zu_res_mean   = 0._wp
1764         zv_res_mean   = 0._wp
1765         zvel_res_mean = 0._wp
1766         zmke          = 0._wp
1767
1768      ELSE
1769
1770         ! * Mean residual (N/m2)
1771         ! Here we take the residual of the linear system (N/m2),
1772         ! We define it as in mitgcm: global area-weighted mean of square-root residual
1773         ! Local residual r = Ax - B expresses to which extent the momentum balance is verified
1774         ! i.e., how close we are to a solution
1775         !!!! UNITS OF RESIDUAL ARE NOT m/s BUT N/m2 (CHECK)
1776         zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp
1777
1778         DO jj = 2, jpj - 1
1779            DO ji = 2, jpi - 1                                     
1780               zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               &
1781                  &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
1782               zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               &
1783                  &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
1784
1785!              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)
1786!              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)
1787   
1788               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
1789               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
1790   
1791            END DO
1792         END DO                 
1793
1794         ! Global ice-concentration, grid-cell-area weighted mean
1795         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_res,  'U', 1., zv_res , 'V', 1. )
1796   
1797         zu_res_mean   = glob_sum( 'ice_rhg_vp', zu_res(:,:) )
1798         zv_res_mean   = glob_sum( 'ice_rhg_vp', zv_res(:,:) )
1799         zvel_res_mean = 0.5_wp * ( zu_res_mean + zv_res_mean )
1800   
1801         ! --- Global mean kinetic energy per unit area (J/m2)
1802         zvel2(:,:) = 0._wp
1803         DO jj = 2, jpj - 1
1804            DO ji = 2, jpi - 1                   
1805               zu     = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point
1806               zv     = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) )
1807               zvel2(ji,jj)  = zu*zu + zv*zv              ! square of ice velocity at T-point 
1808            END DO
1809         END DO
1810         
1811         zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area
1812   
1813      ENDIF ! kitinntot
1814
1815      !--- Spatially-resolved residual at last iteration to send back to main routine (last iteration only)
1816      !--- Calculation @T-point
1817
1818      IF ( kitinntot == kitinntotmax) THEN
1819
1820         DO jj = 2, jpj - 1 ! @ residuals at u and v points
1821            DO ji = 2, jpi - 1
1822
1823               zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)               &
1824                  &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
1825               zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               &
1826                  &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
1827
1828               zu_res(ji,jj)  = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) 
1829               zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) 
1830
1831            END DO
1832         END DO
1833
1834         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_res,  'U',  1., zv_res , 'V',  1. )
1835
1836         DO jj = 2, jpj - 1         ! average @T-point
1837            DO ji = 2, jpi - 1
1838               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) )
1839            END DO
1840         END DO
1841
1842         CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. )
1843
1844      ELSE
1845
1846         pvel_res(:,:) = 0._wp
1847
1848      ENDIF
1849                 
1850      !                                                ! ==================== !
1851
1852      it_inn_file =  ( kt - 1 ) * kitinntotmax + kitinntot ! time step in the file
1853      it_out_file =  ( kt - 1 ) * kitoutmax  + kitout
1854
1855      ! write variables
1856      IF( lwm ) THEN
1857
1858         istatus = NF90_PUT_VAR( ncvgid, nvarid_ures  , (/zu_res_mean/), (/it_inn_file/), (/1/) )        ! Residuals of the linear system, area weighted mean
1859         istatus = NF90_PUT_VAR( ncvgid, nvarid_vres  , (/zv_res_mean/), (/it_inn_file/), (/1/) )        !
1860         istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvel_res_mean/), (/it_inn_file/), (/1/) )      !
1861
1862         istatus = NF90_PUT_VAR( ncvgid, nvarid_uerr_max  , (/puerr_max/), (/it_inn_file/), (/1/) )      ! Max velocit_inn_filey error, sub-it_inn_fileerates
1863         istatus = NF90_PUT_VAR( ncvgid, nvarid_verr_max  , (/pverr_max/), (/it_inn_file/), (/1/) )      !
1864         istatus = NF90_PUT_VAR( ncvgid, nvarid_velerr_max, (/zvel_err_max/), (/it_inn_file/), (/1/) )   !
1865
1866         istatus = NF90_PUT_VAR( ncvgid, nvarid_umad    , (/zu_mad/)  , (/it_inn_file/), (/1/) )         ! velocit_inn_filey MAD, area/sic-weighted, sub-it_inn_fileerates
1867         istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad    , (/zv_mad/)  , (/it_inn_file/), (/1/) )         !
1868         istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad  , (/zvel_mad/), (/it_inn_file/), (/1/) )         !
1869
1870         istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/kitinntot/), (/1/) )                    ! mean kinetic energy
1871
1872         IF ( kitinn == kitinnmax ) THEN ! only print outer mad at the end of inner loop
1873
1874            istatus = NF90_PUT_VAR( ncvgid, nvarid_umad_outer    , (/zu_mad_outer/)  , (/it_out_file/), (/1/) )   ! velocity MAD, area/sic-weighted, outer-iterates
1875            istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad_outer    , (/zv_mad_outer/)  , (/it_out_file/), (/1/) )   !
1876            istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad_outer  , (/zvel_mad_outer/), (/it_out_file/), (/1/) )   !
1877
1878         ENDIF
1879
1880         ! spatially-dependent things (to be coded)
1881         ! velocity difference ?
1882         ! residual ?
1883!        !
1884!        IF ( kiter == kitermax ) THEN
1885!           WRITE(numout,*) ' Should plot the spatially dependent residual '
1886!           istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) )          ! U-residual, spatially dependent
1887!           istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) )          ! V-residual, spatially dependent
1888!        ENDIF
1889
1890         ! close file
1891         IF( kt == nitend )   istatus = NF90_CLOSE( ncvgid )
1892      ENDIF
1893     
1894   END SUBROUTINE rhg_cvg_vp
1895   
1896
1897   
1898#else
1899   !!----------------------------------------------------------------------
1900   !!   Default option         Empty module           NO SI3 sea-ice model
1901   !!----------------------------------------------------------------------
1902#endif
1903
1904   !!==============================================================================
1905END MODULE icedyn_rhg_vp
1906
Note: See TracBrowser for help on using the repository browser.