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

Last change on this file since 13591 was 13591, checked in by vancop, 4 months ago

some debug

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