New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_rhg_vp.F90 in NEMO/trunk/src/ICE – NEMO

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

Last change on this file since 15292 was 15292, checked in by clem, 9 months ago

commit to remove debugging tests in VP rheology. This rheology does not work anyway for the moment

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