source: NEMO/trunk/src/ICE/icedyn_rhg_vp.F90 @ 14433

Last change on this file since 14433 was 14433, checked in by smasson, 3 months ago

trunk: merge dev_r14312_MPI_Interface into the trunk, #2598

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