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

Last change on this file since 15130 was 15014, checked in by smasson, 3 years ago
trunk: simplify F point halo computation, #2693
File size: 85.6 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( 1, 1, 1, 1 )
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( 1, 1, 1, 1 )
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	zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:) ! Snow and ice mass at T-point
334	zmf(:,:) = zmt(:,:) * ff_t(:,:) ! Coriolis factor at T points (m*f)
335
336	DO jj = 2, jpj - 1
337	DO ji = 2, jpi - 1
338
339	! Ice fraction at U-V points
340	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)
341	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)
342
343	! Snow and ice mass at U-V points
344	zm1 = zmt(ji,jj)
345	zm2 = zmt(ji+1,jj)
346	zm3 = zmt(ji,jj+1)
347
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 DO
392	END DO
393
394	CALL iom_put( 'zmsk00x' , zmsk00x ) ! MV DEBUG
395	CALL iom_put( 'zmsk00y' , zmsk00y ) ! MV DEBUG
396	CALL iom_put( 'zmsk01x' , zmsk01x ) ! MV DEBUG
397	CALL iom_put( 'zmsk01y' , zmsk01y ) ! MV DEBUG
398	CALL iom_put( 'ztaux_ai' , ztaux_ai ) ! MV DEBUG
399	CALL iom_put( 'ztauy_ai' , ztauy_ai ) ! MV DEBUG
400	CALL iom_put( 'zspgU' , zspgU ) ! MV DEBUG
401	CALL iom_put( 'zspgV' , zspgV ) ! MV DEBUG
402
403	!------------------------------------------------------------------------------!
404	!
405	! --- Start outer loop
406	!
407	!------------------------------------------------------------------------------!
408
409	zu_c(:,:) = u_ice(:,:)
410	zv_c(:,:) = v_ice(:,:)
411
412	jter = 0
413
414	DO i_out = 1, nn_vp_nout
415
416	IF( lwp ) WRITE(numout,*) ' outer loop i_out : ', i_out
417
418	! Velocities used in the non linear terms are the average of the past two iterates
419	! u_it = 0.5 * ( u_{it-1} + u_{it-2})
420	! Also used in Hibler and Ackley (1983); Zhang and Hibler (1997); Lemieux and Tremblay (2009)
421	zu_c(:,:) = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) )
422	zv_c(:,:) = 0.5_wp * ( v_ice(:,:) + zv_c(:,:) )
423
424	!------------------------------------------------------------------------------!
425	!
426	! --- Right-hand side (RHS) of the linear problem
427	!
428	!------------------------------------------------------------------------------!
429	! In the outer loop, one needs to update all RHS terms
430	! with explicit velocity dependencies (viscosities, coriolis, ocean stress)
431	! as a function of uc
432	!
433
434	!------------------------------------------
435	! -- Strain rates, viscosities and P/Delta
436	!------------------------------------------
437
438	! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
439	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
440	DO ji = 1, jpi - 1
441
442	! shear at F points
443	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) &
444	& + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
445	& ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
446
447	END DO
448	END DO
449
450	CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan
451	CALL iom_put( 'zds' , zds ) ! MV DEBUG
452
453	IF( lwp ) WRITE(numout,*) ' outer loop 1a i_out : ', i_out
454
455	!DO jj = 2, jpj - 1 ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
456	! DO ji = 2, jpi - 1 !
457
458	! MV DEBUG
459	DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
460	DO ji = 2, jpi !
461	! END MV DEBUG
462
463	! shear**2 at T points (doc eq. A16)
464	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
465	& + 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) &
466	& ) * 0.25_wp * r1_e1e2t(ji,jj)
467
468	! divergence at T points
469	zdiv = ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj) &
470	& + e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1) &
471	& ) * r1_e1e2t(ji,jj)
472	zdiv2 = zdiv * zdiv
473
474	! tension at T points
475	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) &
476	& - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
477	& ) * r1_e1e2t(ji,jj)
478	zdt2 = zdt * zdt
479
480	! delta at T points
481	zdeltat = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
482
483	! delta* at T points (following Lemieux and Dupont, GMD 2020)
484	zdeltastar_t(ji,jj) = zdeltat + rn_creepl
485
486	! P/delta at T-points
487	zp_deltastar_t(ji,jj) = strength(ji,jj) / zdeltastar_t(ji,jj)
488
489	! Temporary zzt and zet factors at T-points
490	zzt(ji,jj) = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj)
491	zet(ji,jj) = zzt(ji,jj) * z1_ecc2
492
493	END DO
494	END DO
495
496	CALL lbc_lnk( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. )
497
498	CALL iom_put( 'zzt' , zzt ) ! MV DEBUG
499	CALL iom_put( 'zet' , zet ) ! MV DEBUG
500	CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG
501
502	IF( lwp ) WRITE(numout,*) ' outer loop 1b i_out : ', i_out
503
504	DO jj = 1, jpj - 1
505	DO ji = 1, jpi - 1
506
507	! P/delta* at F points
508	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) )
509
510	! Temporary zef factor at F-point
511	zef(ji,jj) = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * fimask(ji,jj)
512
513	END DO
514	END DO
515
516	CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. )
517	CALL iom_put( 'zef' , zef ) ! MV DEBUG
518	IF( lwp ) WRITE(numout,*) ' outer loop 1c i_out : ', i_out
519
520	!---------------------------------------------------
521	! -- Ocean-ice drag and Coriolis RHS contributions
522	!---------------------------------------------------
523
524	DO jj = 2, jpj - 1
525	DO ji = 2, jpi - 1
526
527	!--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
528	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)
529	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)
530
531	!--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
532	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) ) &
533	& + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
534	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) ) &
535	& + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
536
537	!--- Ocean-ice drag contributions to RHS
538	ztaux_oi_rhsu(ji,jj) = zCwU(ji,jj) * u_oce(ji,jj)
539	ztauy_oi_rhsv(ji,jj) = zCwV(ji,jj) * v_oce(ji,jj)
540
541	! --- U-component of Coriolis Force (energy conserving formulation)
542	! Note Lemieux et al 2008 recommend to do that implicitly, but I don't really see how this could be done
543	zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &
544	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * zv_c(ji ,jj) + e1v(ji ,jj-1) * zv_c(ji ,jj-1) ) &
545	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * zv_c(ji+1,jj) + e1v(ji+1,jj-1) * zv_c(ji+1,jj-1) ) )
546
547	! --- V-component of Coriolis Force (energy conserving formulation)
548	zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &
549	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * zu_c(ji,jj ) + e2u(ji-1,jj ) * zu_c(ji-1,jj ) ) &
550	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) )
551
552	END DO
553	END DO
554
555	IF( lwp ) WRITE(numout,*) ' outer loop 1d i_out : ', i_out
556
557	CALL lbc_lnk( 'icedyn_rhg_vp', zCwU , 'U', -1., zCwV, 'V', -1. )
558	CALL lbc_lnk( 'icedyn_rhg_vp', zCorU, 'U', -1., zCorV, 'V', -1. )
559
560	CALL iom_put( 'zCwU' , zCwU ) ! MV DEBUG
561	CALL iom_put( 'zCwV' , zCwV ) ! MV DEBUG
562	CALL iom_put( 'zCorU' , zCorU ) ! MV DEBUG
563	CALL iom_put( 'zCorV' , zCorV ) ! MV DEBUG
564
565	IF( lwp ) WRITE(numout,*) ' outer loop 1f i_out : ', i_out
566
567	! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system,
568	! so there is no need for lbclnk on drag and coriolis
569
570	!-------------------------------------
571	! -- Internal stress RHS contribution
572	!-------------------------------------
573
574	! --- Stress contributions at T-points
575	DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
576	DO ji = 2, jpi !
577
578	! sig1 contribution to RHS of U-equation at T-points
579	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 )
580
581	! sig2 contribution to RHS of U-equation at T-points
582	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)
583
584	! sig1 contribution to RHS of V-equation at T-points
585	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 )
586
587	! sig2 contribution to RHS of V-equation at T-points
588	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)
589
590	END DO
591	END DO
592
593	CALL iom_put( 'zs1_rhsu' , zs1_rhsu ) ! MV DEBUG
594	CALL iom_put( 'zs2_rhsu' , zs2_rhsu ) ! MV DEBUG
595	CALL iom_put( 'zs1_rhsv' , zs1_rhsv ) ! MV DEBUG
596	CALL iom_put( 'zs2_rhsv' , zs2_rhsv ) ! MV DEBUG
597
598	! a priori, no lbclnk, because rhsu is only used in the inner domain
599
600	! --- Stress contributions at f-points
601	! MV NOTE: I applied fimask on zds, by mimetism on EVP, but without deep understanding of what I was doing
602	! My guess is that this is the way to enforce boundary conditions on strain rate tensor
603
604	IF( lwp ) WRITE(numout,*) ' outer loop 2 i_out : ', i_out
605
606	DO jj = 1, jpj - 1
607	DO ji = 1, jpi - 1
608
609	! sig12 contribution to RHS of U equation at F-points
610	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)
611
612	! sig12 contribution to RHS of V equation at F-points
613	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)
614
615	END DO
616	END DO
617
618	CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. )
619	CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. )
620
621	CALL iom_put( 'zs12_rhsu' , zs12_rhsu ) ! MV DEBUG
622	CALL iom_put( 'zs12_rhsv' , zs12_rhsv ) ! MV DEBUG
623
624	! a priori, no lbclnk, because rhsu are only used in the inner domain
625
626	! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002)
627	! OPT: merge with next loop and use intermediate scalars for zf_rhsu
628
629	DO jj = 2, jpj - 1
630	DO ji = 2, jpi - 1
631	! --- U component of internal force contribution to RHS at U points
632	zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * &
633	( e2u(ji,jj) * ( zs1_rhsu(ji+1,jj) - zs1_rhsu(ji,jj) ) &
634	& + 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) ) &
635	& + 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) ) )
636
637	! --- V component of internal force contribution to RHS at V points
638	zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * &
639	& ( e1v(ji,jj) * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) ) &
640	& + 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) ) &
641	& + 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) ) )
642
643	END DO
644	END DO
645
646	CALL iom_put( 'zf_rhsu' , zf_rhsu ) ! MV DEBUG
647	CALL iom_put( 'zf_rhsv' , zf_rhsv ) ! MV DEBUG
648
649	!---------------------------
650	! -- Sum RHS contributions
651	!---------------------------
652	!
653	! OPT: could use intermediate scalars to reduce memory access
654	DO jj = 2, jpj - 1
655	DO ji = 2, jpi - 1
656
657	! still miss ice ocean stress and acceleration contribution
658	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)
659	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)
660
661	END DO
662	END DO
663
664	CALL lbc_lnk( 'icedyn_rhg_vp', zrhsu, 'U', -1., zrhsv, 'V', -1.)
665	CALL lbc_lnk( 'icedyn_rhg_vp', zmU_t, 'U', -1., zmV_t, 'V', -1.)
666	CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V', -1.)
667
668	CALL iom_put( 'zmU_t' , zmU_t ) ! MV DEBUG
669	CALL iom_put( 'zmV_t' , zmV_t ) ! MV DEBUG
670	CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu ) ! MV DEBUG
671	CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv ) ! MV DEBUG
672	CALL iom_put( 'zrhsu' , zrhsu ) ! MV DEBUG
673	CALL iom_put( 'zrhsv' , zrhsv ) ! MV DEBUG
674
675	! inner domain calculations -> no lbclnk
676
677	IF( lwp ) WRITE(numout,*) ' outer loop 4 i_out : ', i_out
678
679	!---------------------------------------------------------------------------------------!
680	!
681	! --- Linear system matrix
682	!
683	!---------------------------------------------------------------------------------------!
684
685	! Linear system matrix contains all implicit contributions
686	! 1) internal forces, 2) acceleration, 3) ice-ocean drag
687
688	! The linear system equation is written as follows
689	! AU * u_{i-1,j} + BU * u_{i,j} + CU * u_{i+1,j}
690	! = DU * u_{i,j-1} + EU * u_{i,j+1} + RHS (! my convention, not the same as ZH97 )
691
692	! MV Note 1: martin losch applies boundary condition to BU in mitGCM - check whether it is necessary here ?
693	! MV Note 2: "T" factor calculations can be optimized by putting things out of the loop
694	! only zzt and zet are iteration-dependent, other only depend on scale factors
695
696	DO ji = 2, jpi - 1 ! internal domain do loop
697	DO jj = 2, jpj - 1
698
699	!-------------------------------------
700	! -- Internal forces LHS contribution
701	!-------------------------------------
702	!
703	! --- U-component
704	!
705	! "T" factors (intermediate results)
706	!
707	zfac = 0.5_wp * r1_e1e2u(ji,jj)
708	zfac1 = zfac * e2u(ji,jj)
709	zfac2 = zfac * r1_e2u(ji,jj)
710	zfac3 = 2._wp * zfac * r1_e1u(ji,jj)
711
712	zt12U = - zfac1 * zzt(ji+1,jj)
713	zt11U = zfac1 * zzt(ji,jj)
714
715	zt22U = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj)
716	zt21U = zfac2 * zet(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj)
717
718	zt122U = - zfac3 * zef(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj)
719	zt121U = zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)
720
721	!
722	! Linear system coefficients
723	!
724	zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj)
725	zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj)
726	zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj)
727
728	zDU(ji,jj) = zt121U * r1_e1u(ji,jj-1)
729	zEU(ji,jj) = zt122U * r1_e1u(ji,jj+1)
730
731	!
732	! --- V-component
733	!
734	! "T" factors (intermediate results)
735	!
736	zfac = 0.5_wp * r1_e1e2v(ji,jj)
737	zfac1 = zfac * e2v(ji,jj)
738	zfac2 = zfac * r1_e1v(ji,jj)
739	zfac3 = 2._wp * zfac * r1_e2v(ji,jj)
740
741	zt12V = - zfac1 * zzt(ji,jj+1)
742	zt11V = zfac1 * zzt(ji,jj)
743
744	zt22V = zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1)
745	zt21V = - zfac2 * zet(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj)
746
747	zt122V = zfac3 * zef(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj)
748	zt121V = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)
749
750	!
751	! Linear system coefficients
752	!
753	zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1)
754	zBV(ji,jj) = ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj)
755	zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1)
756
757	zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here
758	zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj)
759
760	!-----------------------------------------------------
761	! -- Ocean-ice drag and acceleration LHS contribution
762	!-----------------------------------------------------
763	zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj)
764	zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj)
765
766	END DO
767	END DO
768
769	CALL lbc_lnk( 'icedyn_rhg_vp', zAU , 'U', 1., zAV , 'V', 1. )
770	CALL lbc_lnk( 'icedyn_rhg_vp', zBU , 'U', 1., zBV , 'V', 1. )
771	CALL lbc_lnk( 'icedyn_rhg_vp', zCU , 'U', 1., zCV , 'V', 1. )
772	CALL lbc_lnk( 'icedyn_rhg_vp', zDU , 'U', 1., zDV , 'V', 1. )
773	CALL lbc_lnk( 'icedyn_rhg_vp', zEU , 'U', 1., zEV , 'V', 1. )
774
775	CALL iom_put( 'zAU' , zAU ) ! MV DEBUG
776	CALL iom_put( 'zBU' , zBU ) ! MV DEBUG
777	CALL iom_put( 'zCU' , zCU ) ! MV DEBUG
778	CALL iom_put( 'zDU' , zDU ) ! MV DEBUG
779	CALL iom_put( 'zEU' , zEU ) ! MV DEBUG
780	CALL iom_put( 'zAV' , zAV ) ! MV DEBUG
781	CALL iom_put( 'zBV' , zBV ) ! MV DEBUG
782	CALL iom_put( 'zCV' , zCV ) ! MV DEBUG
783	CALL iom_put( 'zDV' , zDV ) ! MV DEBUG
784	CALL iom_put( 'zEV' , zEV ) ! MV DEBUG
785
786	!------------------------------------------------------------------------------!
787	!
788	! --- Inner loop: solve linear system, check convergence
789	!
790	!------------------------------------------------------------------------------!
791
792	! Inner loop solves the linear problem .. requires 1500 iterations
793	ll_u_iterate = .TRUE.
794	ll_v_iterate = .TRUE.
795
796	DO i_inn = 1, nn_vp_ninn ! inner loop iterations
797
798	IF( lwp ) WRITE(numout,*) ' inner loop 1 i_inn : ', i_inn
799
800	!--- mitgcm computes initial value of residual here...
801
802	jter = jter + 1
803	! l_full_nf_update = jter == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
804
805	zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate
806	zv_b(:,:) = v_ice(:,:)
807
808	! zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp
809	! zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp
810
811	IF ( ll_u_iterate .OR. ll_v_iterate ) THEN
812
813	! ---------------------------- !
814	IF ( ll_u_iterate ) THEN ! --- Solve for u-velocity --- !
815	! ---------------------------- !
816
817	! What follows could be subroutinized...
818
819	! Thomas Algorithm for tridiagonal solver
820	! Au(i-1,j)+Bu(i,j)+C*u(i+1,j) = F
821
822	zFU(:,:) = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; zCU_prime(:,:) = 0._wp
823
824	DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! )
825
826	! OPT: could be even better optimized with a true red-black SOR
827
828	IF ( jn == 1 ) THEN ; jj_min = 2
829	ELSE ; jj_min = 3
830	ENDIF
831
832	IF ( lwp ) WRITE(numout,*) ' Into the U-zebra loop at step jn = ', jn, ', with jj_min = ', jj_min
833
834	DO jj = jj_min, jpj - 1, nn_zebra_vp
835
836	!------------------------
837	! Independent term (zFU)
838	!------------------------
839	DO ji = 2, jpi - 1
840
841	! boundary condition substitution
842	! see Zhang and Hibler, 1997, Appendix B
843	zAA3 = 0._wp
844	IF ( ji == 2 ) zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj)
845	IF ( ji == jpi - 1 ) zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj)
846
847	! right hand side
848	zFU(ji,jj) = ( zrhsu(ji,jj) & ! right-hand side terms
849	& + zAA3 & ! boundary condition translation
850	& + zDU(ji,jj) * u_ice(ji,jj-1) & ! internal force, j-1
851	& + zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1
852
853	END DO
854
855	END DO
856
857	CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', 1. )
858
859	!---------------
860	! Forward sweep
861	!---------------
862	DO jj = jj_min, jpj - 1, nn_zebra_vp
863
864	DO ji = 3, jpi - 1
865
866	zfac = SIGN( 1._wp , zBU(ji-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU(ji-1,jj) ) - epsi20 ) )
867	zw = zfac * zAU(ji,jj) / MAX ( ABS( zBU(ji-1,jj) ) , epsi20 )
868	zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj)
869	zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj)
870
871	END DO
872
873	END DO
874
875	CALL lbc_lnk( 'icedyn_rhg_vp', zFU_prime, 'U', 1., zBU_prime, 'U', 1. )
876
877	!-----------------------------
878	! Backward sweep & relaxation
879	!-----------------------------
880
881	DO jj = jj_min, jpj - 1, nn_zebra_vp
882
883	! --- Backward sweep
884	! last row
885	zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) )
886	u_ice(jpi-1,jj) = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) &
887	& * umask(jpi-1,jj,1)
888	DO ji = jpi-2 , 2, -1 ! all other rows ! ---> original backward loop
889	zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) )
890	u_ice(ji,jj) = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1) &
891	& / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 )
892	END DO
893
894	!--- Relaxation
895	! and velocity masking for little-ice and no-ice cases
896	DO ji = 2, jpi - 1
897
898	u_ice(ji,jj) = zu_b(ji,jj) + zrelaxu_vp * ( u_ice(ji,jj) - zu_b(ji,jj) ) ! relaxation
899
900	u_ice(ji,jj) = zmsk00x(ji,jj) & ! masking
901	& * ( zmsk01x(ji,jj) * u_ice(ji,jj) &
902	& + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp ) * umask(ji,jj,1)
903
904	END DO
905
906	END DO ! jj
907
908	END DO ! zebra loop
909
910	ENDIF ! ll_u_iterate
911
912	! ! ---------------------------- !
913	IF ( ll_v_iterate ) THEN ! --- Solve for V-velocity --- !
914	! ! ---------------------------- !
915
916	! MV OPT: what follows could be subroutinized...
917	! Thomas Algorithm for tridiagonal solver
918	! Av(i,j-1)+Bv(i,j)+C*v(i,j+1) = F
919	! It is intentional to have a ji then jj loop for V-velocity
920	!!! ZH97 explain it is critical for convergence speed
921
922	zFV(:,:) = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; zCV_prime(:,:) = 0._wp
923
924	DO jn = 1, nn_zebra_vp ! "zebra" loop
925
926	IF ( jn == 1 ) THEN ; ji_min = 2
927	ELSE ; ji_min = 3
928	ENDIF
929
930	IF ( lwp ) WRITE(numout,*) ' Into the V-zebra loop at step jn = ', jn, ', with ji_min = ', ji_min
931
932	DO ji = ji_min, jpi - 1, nn_zebra_vp
933
934	!------------------------
935	! Independent term (zFV)
936	!------------------------
937	DO jj = 2, jpj - 1
938
939	! boundary condition substitution (check it is correctly applied !!!)
940	! see Zhang and Hibler, 1997, Appendix B
941	zAA3 = 0._wp
942	IF ( jj == 2 ) zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1)
943	IF ( jj == jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1)
944
945	! right hand side
946	zFV(ji,jj) = ( zrhsv(ji,jj) & ! right-hand side terms
947	& + zAA3 & ! boundary condition translation
948	& + zDV(ji,jj) * v_ice(ji-1,jj) & ! internal force, j-1
949	& + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1) ! internal force, j+1
950
951	END DO
952
953	END DO
954
955	CALL lbc_lnk( 'icedyn_rhg_vp', zFV, 'V', 1.)
956
957	!---------------
958	! Forward sweep
959	!---------------
960	DO ji = ji_min, jpi - 1, nn_zebra_vp
961
962	DO jj = 3, jpj - 1
963
964	zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) )
965	zw = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 )
966	zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1)
967	zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1)
968
969	END DO
970
971	END DO
972
973	CALL lbc_lnk( 'icedyn_rhg_vp', zFV_prime, 'V', 1., zBV_prime, 'V', 1. )
974
975	!-----------------------------
976	! Backward sweep & relaxation
977	!-----------------------------
978	DO ji = ji_min, jpi - 1, nn_zebra_vp
979
980	! --- Backward sweep
981	! last row
982	zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) )
983	v_ice(ji,jpj-1) = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) &
984	& * vmask(ji,jpj-1,1) ! last row
985
986	! other rows
987	DO jj = jpj-2, 2, -1 ! original back loop
988	zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) )
989	v_ice(ji,jj) = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) &
990	& / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 )
991	END DO
992
993	! --- Relaxation & masking (should it be now or later)
994	DO jj = 2, jpj - 1
995
996	v_ice(ji,jj) = zv_b(ji,jj) + zrelaxv_vp * ( v_ice(ji,jj) - zv_b(ji,jj) ) ! relaxation
997
998	v_ice(ji,jj) = zmsk00y(ji,jj) & ! masking
999	& * ( zmsk01y(ji,jj) * v_ice(ji,jj) &
1000	& + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp ) * vmask(ji,jj,1)
1001
1002	END DO ! jj
1003
1004	END DO ! ji
1005
1006	END DO ! zebra loop
1007
1008	ENDIF ! ll_v_iterate
1009
1010	CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
1011
1012	!--------------------------------------------------------------------------------------
1013	! -- Check convergence based on maximum velocity difference, continue or stop the loop
1014	!--------------------------------------------------------------------------------------
1015
1016	!------
1017	! on U
1018	!------
1019	! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check
1020	! then we must optimize the use of the mpp_max, which is prohibitive
1021	zuerr_max = 0._wp
1022
1023	IF ( ll_u_iterate .AND. MOD ( i_inn, nn_vp_chkcvg ) == 0 ) THEN
1024
1025	! - Maximum U-velocity difference
1026	zuerr(:,:) = 0._wp
1027	DO jj = 2, jpj - 1
1028	DO ji = 2, jpi - 1
1029	zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1)
1030	END DO
1031	END DO
1032	zuerr_max = MAXVAL( zuerr )
1033	CALL mpp_max( 'icedyn_rhg_evp', zuerr_max ) ! max over the global domain - damned!
1034
1035	! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)
1036	IF ( i_inn > 1 .AND. zuerr_max > zuerr_max_vp ) THEN
1037	IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
1038	ll_u_iterate = .FALSE.
1039	ENDIF
1040
1041	! - Stop if error small enough
1042	IF ( zuerr_max < zuerr_min_vp ) THEN
1043	IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer U-iteration ", i_out, " after ", i_inn, " iterations, finished! "
1044	ll_u_iterate = .FALSE.
1045	ENDIF
1046
1047	ENDIF ! ll_u_iterate
1048
1049	!------
1050	! on V
1051	!------
1052	zverr_max = 0._wp
1053
1054	IF ( ll_v_iterate .AND. MOD ( i_inn, nn_vp_chkcvg ) == 0 ) THEN
1055
1056	! - Maximum V-velocity difference
1057	zverr(:,:) = 0._wp
1058	DO jj = 2, jpj - 1
1059	DO ji = 2, jpi - 1
1060	zverr(ji,jj) = ABS ( ( v_ice(ji,jj) - zv_b(ji,jj) ) ) * vmask(ji,jj,1)
1061	END DO
1062	END DO
1063
1064	zverr_max = MAXVAL( zverr )
1065	CALL mpp_max( 'icedyn_rhg_evp', zverr_max ) ! max over the global domain - damned!
1066
1067	! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)
1068	IF ( i_inn > 1 .AND. zverr_max > zuerr_max_vp ) THEN
1069	IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zverr_max, " in outer V-iteration ", i_out, " after ", i_inn, " iterations, we stopped "
1070	ll_v_iterate = .FALSE.
1071	ENDIF
1072
1073	! - Stop if error small enough
1074	IF ( zverr_max < zuerr_min_vp ) THEN
1075	IF ( lwp ) WRITE(numout,*) " VP rheology nicely done in outer V-iteration ", i_out, " after ", i_inn, " iterations, finished! "
1076	ll_v_iterate = .FALSE.
1077	ENDIF
1078
1079	ENDIF ! ll_v_iterate
1080
1081	ENDIF ! --- end ll_u_iterate or ll_v_iterate
1082
1083	!---------------------------------------------------------------------------------------
1084	!
1085	! --- Calculate extra convergence diagnostics and save them
1086	!
1087	!---------------------------------------------------------------------------------------
1088
1089	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, &
1090	& zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV )
1091
1092	IF ( lwp ) WRITE(numout,*) ' Done convergence tests '
1093
1094	END DO ! i_inn, end of inner loop
1095
1096	END DO ! End of outer loop (i_out) =============================================================================================
1097
1098	IF ( lwp ) WRITE(numout,*) ' We are out of outer loop '
1099
1100	CALL lbc_lnk( 'icedyn_rhg_vp', zFU , 'U', 1., zFV , 'V', 1. )
1101	CALL lbc_lnk( 'icedyn_rhg_vp', zBU_prime , 'U', 1., zBV_prime , 'V', 1. )
1102	CALL lbc_lnk( 'icedyn_rhg_vp', zFU_prime , 'U', 1., zFV_prime , 'V', 1. )
1103	CALL lbc_lnk( 'icedyn_rhg_vp', zCU_prime , 'U', 1., zCV_prime , 'V', 1. )
1104
1105	CALL iom_put( 'zFU' , zFU ) ! MV DEBUG
1106	CALL iom_put( 'zBU_prime' , zBU_prime ) ! MV DEBUG
1107	CALL iom_put( 'zCU_prime' , zCU_prime ) ! MV DEBUG
1108	CALL iom_put( 'zFU_prime' , zFU_prime ) ! MV DEBUG
1109
1110	CALL iom_put( 'zFV' , zFV ) ! MV DEBUG
1111	CALL iom_put( 'zBV_prime' , zBV_prime ) ! MV DEBUG
1112	CALL iom_put( 'zCV_prime' , zCV_prime ) ! MV DEBUG
1113	CALL iom_put( 'zFV_prime' , zFV_prime ) ! MV DEBUG
1114
1115	CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
1116
1117	IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg '
1118	IF( iom_use('uice_dbg' ) ) CALL iom_put( 'uice_dbg' , u_ice ) ! ice velocity u after solver
1119	IF( iom_use('vice_dbg' ) ) CALL iom_put( 'vice_dbg' , v_ice ) ! ice velocity v after solver
1120
1121	!------------------------------------------------------------------------------!
1122	!
1123	! --- Convergence diagnostics
1124	!
1125	!------------------------------------------------------------------------------!
1126
1127	IF( nn_rhg_chkcvg /= 0 ) THEN
1128
1129	IF( iom_use('uice_cvg') ) THEN
1130	CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , & ! ice velocity difference at last iteration
1131	& ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
1132	ENDIF
1133
1134	ENDIF
1135
1136	! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead
1137	DO jj = 2, jpj - 1
1138	DO ji = 2, jpi - 1
1139
1140	u_ice(ji,jj) = zmsk00x(ji,jj) &
1141	& * ( zmsk01x(ji,jj) * u_oce(ji,jj) * 0.01_wp &
1142	+ ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp )
1143
1144	v_ice(ji,jj) = zmsk00y(ji,jj) &
1145	& * ( zmsk01y(ji,jj) * v_oce(ji,jj) * 0.01_wp &
1146	+ ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp )
1147
1148	END DO
1149	END DO
1150
1151	CALL lbc_lnk( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )
1152
1153	IF ( lwp ) WRITE(numout,*) ' Velocity replaced '
1154
1155	! END DEBUG
1156
1157	!------------------------------------------------------------------------------!
1158	!
1159	! --- Recompute delta, shear and div (inputs for mechanical redistribution)
1160	!
1161	!------------------------------------------------------------------------------!
1162	!
1163	! MV OPT: subroutinize ?
1164
1165	DO jj = 1, jpj - 1
1166	DO ji = 1, jpi - 1
1167
1168	! shear at F points
1169	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) &
1170	& + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) &
1171	& ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
1172
1173	END DO
1174	END DO
1175
1176	DO jj = 2, jpj - 1
1177	DO ji = 2, jpi - 1 !
1178
1179	! tension**2 at T points
1180	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) &
1181	& - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &
1182	& ) * r1_e1e2t(ji,jj)
1183	zdt2 = zdt * zdt
1184
1185	zten_i(ji,jj) = zdt
1186
1187	! shear**2 at T points (doc eq. A16)
1188	zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &
1189	& + 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) &
1190	& ) * 0.25_wp * r1_e1e2t(ji,jj)
1191
1192	! shear at T points
1193	pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
1194
1195	! divergence at T points
1196	pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) &
1197	& + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) &
1198	& ) * r1_e1e2t(ji,jj)
1199
1200	! delta at T points
1201	zdelta = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )
1202	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
1203
1204	!pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
1205	pdelta_i(ji,jj) = zdelta + rn_creepl
1206
1207	END DO
1208	END DO
1209
1210	IF ( lwp ) WRITE(numout,*) ' Deformation recalculated '
1211
1212	CALL lbc_lnk( 'icedyn_rhg_vp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
1213
1214	!------------------------------------------------------------------------------!
1215	!
1216	! --- Diagnostics
1217	!
1218	!------------------------------------------------------------------------------!
1219	!
1220	! MV OPT: subroutinize ?
1221	!
1222	!----------------------------------
1223	! --- Recompute stresses if needed
1224	!----------------------------------
1225	!
1226	! ---- Sea ice stresses at T-points
1227	IF ( iom_use('normstr') .OR. iom_use('sheastr') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1228
1229	DO jj = 2, jpj - 1
1230	DO ji = 2, jpi - 1
1231	zp_deltastar_t(ji,jj) = strength(ji,jj) / pdelta_i(ji,jj)
1232	zfac = zp_deltastar_t(ji,jj)
1233	zs1(ji,jj) = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
1234	zs2(ji,jj) = zfac * z1_ecc2 * zten_i(ji,jj)
1235	zs12(ji,jj) = zfac * z1_ecc2 * pshear_i(ji,jj)
1236	END DO
1237	END DO
1238
1239	CALL lbc_lnk( 'icedyn_rhg_vp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'T', 1. )
1240
1241	ENDIF
1242
1243	! ---- s12 at F-points
1244	IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1245
1246	DO jj = 1, jpj - 1
1247	DO ji = 1, jpi - 1
1248
1249	! P/delta* at F points
1250	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) )
1251
1252	! s12 at F-points
1253	zs12f(ji,jj) = zp_deltastar_f * z1_ecc2 * zds(ji,jj)
1254
1255	END DO
1256	END DO
1257
1258	CALL lbc_lnk( 'icedyn_rhg_vp', zs12f, 'F', 1. )
1259
1260	ENDIF
1261
1262	IF ( lwp ) WRITE(numout,*) ' zs12f recalculated '
1263
1264	!
1265	!-----------------------
1266	! --- Store diagnostics
1267	!-----------------------
1268	!
1269	! --- Ice-ocean, ice-atm. & ice-ocean bottom (landfast) stresses --- !
1270	IF( iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
1271	& iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
1272
1273	ALLOCATE( ztaux_oi(jpi,jpj) , ztauy_oi(jpi,jpj) )
1274
1275	!--- Recalculate oceanic stress at last inner iteration
1276	DO jj = 2, jpj - 1
1277	DO ji = 2, jpi - 1
1278
1279	!--- ice u-velocity @V points, v-velocity @U points (for non-linear drag computation)
1280	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)
1281	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)
1282
1283	!--- non-linear drag coefficients (need to be updated at each outer loop, see Lemieux and Tremblay JGR09, p.3, beginning of Section 3)
1284	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) ) &
1285	& + ( zv_cU - v_oceU(ji,jj) ) * ( zv_cU - v_oceU(ji,jj) ) )
1286	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) ) &
1287	& + ( zu_cV - u_oceV(ji,jj) ) * ( zu_cV - u_oceV(ji,jj) ) )
1288
1289	!--- Ocean-ice stress
1290	ztaux_oi(ji,jj) = zCwU(ji,jj) * ( u_oce(ji,jj) - u_ice(ji,jj) )
1291	ztauy_oi(ji,jj) = zCwV(ji,jj) * ( v_oce(ji,jj) - v_ice(ji,jj) )
1292
1293	END DO
1294	END DO
1295
1296	!
1297	CALL lbc_lnk( 'icedyn_rhg_vp', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) !, &
1298	! & ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
1299	!
1300	CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
1301	CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
1302	CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
1303	CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
1304	! CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
1305	! CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
1306
1307	DEALLOCATE( ztaux_oi , ztauy_oi )
1308
1309	ENDIF
1310
1311	! --- Divergence, shear and strength --- !
1312	IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence
1313	IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear
1314	IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta
1315	IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength
1316
1317	IF ( lwp ) WRITE(numout,*) 'Some terms recalculated '
1318
1319	! --- Stress tensor invariants (SIMIP diags) --- !
1320	IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
1321	!
1322	! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags.
1323	! Definitions following Coon (1974) and Feltham (2008)
1324	!
1325	! sigma1, sigma2, sigma12 are useful (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
1326	! however these are NOT stress tensor components, neither stress invariants, nor stress principal components
1327	!
1328	ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
1329	!
1330	DO jj = 2, jpj - 1
1331	DO ji = 2, jpi - 1
1332	! Stress invariants
1333	zsig_I(ji,jj) = zs1(ji,jj) * 0.5_wp ! 1st invariant, aka average normal stress aka negative pressure
1334	zsig_II(ji,jj) = SQRT ( zs2(ji,jj) * zs2(ji,jj) * 0.25_wp + zs12(ji,jj) ) ! 2nd invariant, aka maximum shear stress
1335	END DO
1336	END DO
1337
1338	CALL lbc_lnk( 'icedyn_rhg_vp', zsig_I, 'T', 1., zsig_II, 'T', 1.)
1339
1340	IF( iom_use('normstr') ) CALL iom_put( 'normstr' , zsig_I(:,:) * zmsk00(:,:) ) ! Normal stress
1341	IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
1342
1343	DEALLOCATE ( zsig_I, zsig_II )
1344
1345	ENDIF
1346
1347	IF ( lwp ) WRITE(numout,*) 'SIMIP work done'
1348
1349	! --- Normalized stress tensor principal components --- !
1350	! These are used to plot the normalized yield curve (Lemieux & Dupont, GMD 2020)
1351	! To plot the yield curve and evaluate physical convergence, they have two recommendations
1352	! Recommendation 1 : Use ice strength, not replacement pressure
1353	! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765)
1354	! R2 means we need to recompute stresses
1355
1356	IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
1357	!
1358	ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
1359	!
1360	DO jj = 2, jpj - 1
1361	DO ji = 2, jpi - 1
1362	! Ice stresses computed with viscosities (delta, p/delta) at previous iterates
1363	! and deformations at current iterates
1364	! following Lemieux & Dupont (2020)
1365	zfac = zp_deltastar_t(ji,jj)
1366	zsig1 = zfac * ( pdivu_i(ji,jj) - zdeltastar_t(ji,jj) )
1367	zsig1 = 0._wp !!! FUCKING DEBUG TEST !!!
1368	zsig2 = zfac * z1_ecc2 * zten_i(ji,jj)
1369	zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)
1370
1371	! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
1372	zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st invariant
1373	zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 ) ! 2nd invariant
1374
1375	! Normalized principal stresses (used to display the ellipse)
1376	z1_strength = 1._wp / MAX ( 1._wp , strength(ji,jj) )
1377	zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
1378	zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
1379	END DO
1380	END DO
1381	IF ( lwp ) WRITE(numout,*) 'Some shitty stress work done'
1382	!
1383	CALL lbc_lnk( 'icedyn_rhg_vp', zsig1_p, 'T', 1., zsig2_p, 'T', 1.)
1384	!
1385	IF ( lwp ) WRITE(numout,*) ' Beauaaaarflblbllll '
1386	!
1387	CALL iom_put( 'sig1_pnorm' , zsig1_p )
1388	CALL iom_put( 'sig2_pnorm' , zsig2_p )
1389
1390	DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II )
1391
1392	IF ( lwp ) WRITE(numout,*) ' So what ??? '
1393
1394	ENDIF
1395
1396	! --- SIMIP, terms of tendency for momentum equation --- !
1397	IF( iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
1398	& iom_use('corstrx') .OR. iom_use('corstry') ) THEN
1399
1400	! --- Recalculate Coriolis stress at last inner iteration
1401	DO jj = 2, jpj - 1
1402	DO ji = 2, jpi - 1
1403	! --- U-component
1404	zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * &
1405	& ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) &
1406	& + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
1407	zCorV(ji,jj) = - 0.25_wp * r1_e2v(ji,jj) * &
1408	& ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) &
1409	& + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
1410	END DO
1411	END DO
1412	!
1413	CALL lbc_lnk( 'icedyn_rhg_vp', zspgU, 'U', -1., zspgV, 'V', -1., &
1414	& zCorU, 'U', -1., zCorV, 'V', -1. )
1415	!
1416	CALL iom_put( 'dssh_dx' , zspgU * zmsk00 ) ! Sea-surface tilt term in force balance (x)
1417	CALL iom_put( 'dssh_dy' , zspgV * zmsk00 ) ! Sea-surface tilt term in force balance (y)
1418	CALL iom_put( 'corstrx' , zCorU * zmsk00 ) ! Coriolis force term in force balance (x)
1419	CALL iom_put( 'corstry' , zCorV * zmsk00 ) ! Coriolis force term in force balance (y)
1420
1421	ENDIF
1422
1423	IF ( iom_use('intstrx') .OR. iom_use('intstry') ) THEN
1424
1425	! Recalculate internal forces (divergence of stress tensor) at last inner iteration
1426	DO jj = 2, jpj - 1
1427	DO ji = 2, jpi - 1
1428	zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) &
1429	& + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) &
1430	& ) * r1_e2u(ji,jj) &
1431	& + ( zs12f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) &
1432	& ) * 2._wp * r1_e1u(ji,jj) &
1433	& ) * r1_e1e2u(ji,jj)
1434	zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) &
1435	& - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) &
1436	& ) * r1_e1v(ji,jj) &
1437	& + ( zs12f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) &
1438	& ) * 2._wp * r1_e2v(ji,jj) &
1439	& ) * r1_e1e2v(ji,jj)
1440	END DO
1441	END DO
1442
1443	CALL lbc_lnk( 'icedyn_rhg_vp', zfU, 'U', -1., zfV, 'V', -1. )
1444
1445	CALL iom_put( 'intstrx' , zfU * zmsk00 ) ! Internal force term in force balance (x)
1446	CALL iom_put( 'intstry' , zfV * zmsk00 ) ! Internal force term in force balance (y)
1447
1448	ENDIF
1449
1450	! --- Ice & snow mass and ice area transports
1451	IF( iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
1452	& iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
1453	!
1454	ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
1455	& zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
1456	!
1457	DO jj = 2, jpj - 1
1458	DO ji = 2, jpi - 1
1459	! 2D ice mass, snow mass, area transport arrays (X, Y)
1460	zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
1461	zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
1462
1463	zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
1464	zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) ! '' Y- ''
1465
1466	zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
1467	zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) ! '' Y- ''
1468
1469	zdiag_xatrp(ji,jj) = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport, X-component
1470	zdiag_yatrp(ji,jj) = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) ) ! '' Y- ''
1471
1472	END DO
1473	END DO
1474
1475	CALL lbc_lnk( 'icedyn_rhg_vp', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
1476	& zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
1477	& zdiag_xatrp , 'U', -1., zdiag_yatrp , 'V', -1. )
1478
1479	CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice ) ! X-component of sea-ice mass transport (kg/s)
1480	CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice ) ! Y-component of sea-ice mass transport
1481	CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw ) ! X-component of snow mass transport (kg/s)
1482	CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw ) ! Y-component of snow mass transport
1483	CALL iom_put( 'xatrp' , zdiag_xatrp ) ! X-component of ice area transport
1484	CALL iom_put( 'yatrp' , zdiag_yatrp ) ! Y-component of ice area transport
1485
1486	DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
1487	& zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
1488
1489	ENDIF
1490
1491	END SUBROUTINE ice_dyn_rhg_vp
1492
1493
1494
1495	SUBROUTINE rhg_cvg_vp( kt, kiter, kitermax, pu, pv, pmt, puerr_max, pverr_max, pglob_area, &
1496	& prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV )
1497
1498	!!----------------------------------------------------------------------
1499	!! * ROUTINE rhg_cvg_vp *
1500	!!
1501	!! ** Purpose : check convergence of VP ice rheology
1502	!!
1503	!! ** Method : create a file ice_cvg.nc containing a few convergence diagnostics
1504	!! This routine is called every sub-iteration, so it is cpu expensive
1505	!!
1506	!! Calculates / stores
1507	!! - maximum absolute U-V difference (uice_cvg, u_dif, v_dif, m/s)
1508	!! - 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)
1509	!! - mean kinetic energy (mke_ice, J/m2)
1510	!!
1511	!! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)
1512	!!----------------------------------------------------------------------
1513	INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index
1514	REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pmt ! now velocity and mass per unit area
1515	REAL(wp), INTENT(in) :: puerr_max, pverr_max ! absolute mean velocity difference
1516	REAL(wp), INTENT(in) :: pglob_area ! global ice area
1517	REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsu, pAU, pBU, pCU, pDU, pEU ! linear system coefficients
1518	REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsv, pAV, pBV, pCV, pDV, pEV
1519	!!
1520	INTEGER :: it, idtime, istatus, ix_dim, iy_dim
1521	INTEGER :: ji, jj ! dummy loop indices
1522	REAL(wp) :: zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars
1523	REAL(wp) :: z1_pglob_area
1524	REAL(wp), DIMENSION(jpi,jpj) :: zu_res, zv_res, zvel2 ! local arrays
1525
1526	CHARACTER(len=20) :: clname
1527	!!----------------------------------------------------------------------
1528
1529	IF( lwp ) THEN
1530	WRITE(numout,*)
1531	WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control'
1532	WRITE(numout,*) '~~~~~~~~~~~'
1533	WRITE(numout,*) ' kiter = : ', kiter
1534	WRITE(numout,*) ' kitermax = : ', kitermax
1535	ENDIF
1536
1537	! create file
1538	IF( kt == nit000 .AND. kiter == 1 ) THEN
1539	!
1540	IF( lwm ) THEN
1541
1542	clname = 'ice_cvg.nc'
1543	IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
1544	istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
1545
1546	istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime )
1547	istatus = NF90_DEF_DIM( ncvgid, 'x' , jpi, ix_dim )
1548	istatus = NF90_DEF_DIM( ncvgid, 'y' , jpj, iy_dim )
1549
1550	! i suggest vel_dif instead
1551	istatus = NF90_DEF_VAR( ncvgid, 'u_res' , NF90_DOUBLE , (/ idtime /), nvarid_ures )
1552	istatus = NF90_DEF_VAR( ncvgid, 'v_res' , NF90_DOUBLE , (/ idtime /), nvarid_vres )
1553	istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE , (/ idtime /), nvarid_velres )
1554	istatus = NF90_DEF_VAR( ncvgid, 'u_dif' , NF90_DOUBLE , (/ idtime /), nvarid_udif )
1555	istatus = NF90_DEF_VAR( ncvgid, 'v_dif' , NF90_DOUBLE , (/ idtime /), nvarid_vdif )
1556	istatus = NF90_DEF_VAR( ncvgid, 'vel_dif', NF90_DOUBLE , (/ idtime /), nvarid_veldif )
1557	istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE , (/ idtime /), nvarid_mke )
1558
1559	istatus = NF90_DEF_VAR( ncvgid, 'u_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_ures_xy)
1560	istatus = NF90_DEF_VAR( ncvgid, 'v_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_vres_xy)
1561
1562	istatus = NF90_ENDDEF(ncvgid)
1563
1564	ENDIF
1565	!
1566	ENDIF
1567
1568	IF ( lwp ) WRITE(numout,*) ' File created '
1569
1570	! --- Max absolute velocity difference with previous iterate (zveldif)
1571	zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code
1572	! if puerrmask and pverrmax are masked at 15% (TEST)
1573
1574	! --- Mean residual and kinetic energy
1575	IF ( kiter == 1 ) THEN
1576
1577	zu_res_mean = 0._wp
1578	zv_res_mean = 0._wp
1579	zvelres = 0._wp
1580	zmke = 0._wp
1581
1582	ELSE
1583
1584	! -- Mean residual (N/m^2), zu_res_mean
1585	! Here we take the residual of the linear system (N/m^2),
1586	! We define it as in mitgcm: square-root of area-weighted mean square residual
1587	! Local residual r = Ax - B expresses to which extent the momentum balance is verified
1588	! i.e., how close we are to a solution
1589
1590	IF ( lwp ) WRITE(numout,*) ' TEST 1 '
1591
1592	z1_pglob_area = 1._wp / pglob_area
1593
1594	zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp
1595
1596	DO jj = 2, jpj - 1
1597	DO ji = 2, jpi - 1
1598	zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) &
1599	& - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) )
1600
1601	zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) &
1602	& - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) )
1603
1604	zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * e1e2u(ji,jj) * z1_pglob_area
1605	zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area
1606
1607	END DO
1608	END DO
1609
1610	IF ( lwp ) WRITE(numout,*) ' TEST 2 '
1611	zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) )
1612	zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) )
1613	IF ( lwp ) WRITE(numout,*) ' TEST 3 '
1614	zvelres = 0.5_wp * ( zu_res_mean + zv_res_mean )
1615
1616	IF ( lwp ) WRITE(numout,*) ' TEST 4 '
1617
1618	! -- Global mean kinetic energy per unit area (J/m2)
1619	zvel2(:,:) = 0._wp
1620	DO jj = 2, jpj - 1
1621	DO ji = 2, jpi - 1
1622	zu = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point
1623	zv = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) )
1624	zvel2(ji,jj) = zuzu + zvzv ! square of ice velocity at T-point
1625	END DO
1626	END DO
1627
1628	IF ( lwp ) WRITE(numout,*) ' TEST 5 '
1629
1630	zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area
1631
1632	IF ( lwp ) WRITE(numout,*) ' TEST 6 '
1633
1634	ENDIF ! kiter
1635
1636	! ! ==================== !
1637
1638	! time
1639	it = ( kt - 1 ) * kitermax + kiter
1640
1641
1642	IF( lwm ) THEN
1643	! write variables
1644	istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system
1645	istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system
1646	istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) ! average of u- and v- residuals
1647	istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) ! max U velocity difference, inner iterations
1648	istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) ! max V velocity difference, inner iterations
1649	istatus = NF90_PUT_VAR( ncvgid, nvarid_veldif, (/zveldif/), (/it/), (/1/) ) ! max U or V velocity diff between subiterations
1650	istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) ! mean kinetic energy
1651
1652	!
1653	IF ( kiter == kitermax ) THEN
1654	WRITE(numout,*) ' Should plot the spatially dependent residual '
1655	istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) ) ! U-residual, spatially dependent
1656	istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) ) ! V-residual, spatially dependent
1657	ENDIF
1658
1659	! close file
1660	IF( kt == nitend ) istatus = NF90_CLOSE( ncvgid )
1661	ENDIF
1662
1663	END SUBROUTINE rhg_cvg_vp
1664
1665
1666
1667	#else
1668	!!----------------------------------------------------------------------
1669	!! Default option Empty module NO SI3 sea-ice model
1670	!!----------------------------------------------------------------------
1671	#endif
1672
1673	!!==============================================================================
1674	END MODULE icedyn_rhg_vp
1675

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

Download in other formats: