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

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

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

Download in other formats: