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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_rhg_vp.F90 in NEMO/branches/2020/SI3-03_VP_rheology/src/ICE – NEMO

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

Last change on this file since 15403 was 15403, checked in by vancop, 9 months ago

VP with much less bugs running at least 8 time steps on SAS

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