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

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

Last change on this file since 13536 was 13536, checked in by vancop, 20 months ago

VP rheology, first code version finalised

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