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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icedyn_rhg_vp.F90 in NEMO/branches/2020/SI3-03_VP_rheology/src/ICE – NEMO

Context Navigation

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

Last change on this file since 13544 was 13544, checked in by vancop, 4 years ago
now VP rheology compiles
File size: 73.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: