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

Last change on this file since 13593 was 13593, checked in by vancop, 4 months ago

Modified tridiagonal system

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