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

Last change on this file since 15523 was 15523, checked in by vancop, 8 months ago

Some cleaning of lbclnks in VP routine

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