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

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

icedyn_rhg_vp.F90 in NEMO/trunk/src/ICE – NEMO

Context Navigation

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

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

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

Download in other formats: