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

Last change on this file since 13522 was 13522, checked in by vancop, 4 years ago
Close to compiling VP
File size: 70.0 KB

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

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

Download in other formats: