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

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

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

Context Navigation

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

Last change on this file since 13681 was 13681, checked in by vancop, 3 years ago
Post-gurvanian subroutine volume 2 including several bugfixes for Clement check
File size: 85.8 KB

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

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

Download in other formats: