source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90 @ 8486

Last change on this file since 8486 was 8486, checked in by clem, 3 years ago

changes in style - part1 - (now the code looks better txs to Gurvan's comments)

File size: 42.3 KB
Line 
1MODULE icerhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icerhg_evp  ***
4   !!   Ice rheology : sea ice rheology
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
8   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy
9   !!            3.3  !  2009-05  (G.Garric) addition of the evp cas
10   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation
11   !!            3.5  !  2012-08  (R. Benshila)  AGRIF
12   !!            3.6  !  2016-06  (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
13   !!----------------------------------------------------------------------
14#if defined key_lim3
15   !!----------------------------------------------------------------------
16   !!   'key_lim3'                                      LIM-3 sea-ice model
17   !!----------------------------------------------------------------------
18   !!   ice_rhg_evp       : computes ice velocities
19   !!----------------------------------------------------------------------
20   USE phycst         ! Physical constant
21   USE par_oce        ! Ocean parameters
22   USE dom_oce        ! Ocean domain
23   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
24   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
25   USE ice            ! ice variables
26   USE icerdgrft      ! ice strength
27   !
28   USE lbclnk         ! Lateral Boundary Condition / MPP link
29   USE lib_mpp        ! MPP library
30   USE in_out_manager ! I/O manager
31   USE prtctl         ! Print control
32   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
33#if defined key_agrif
34   USE agrif_lim3_interp
35#endif
36   USE bdy_oce   , ONLY: ln_bdy 
37   USE bdyice 
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   ice_rhg_evp        ! routine called by icerhg.F90
43
44   !! * Substitutions
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
48   !! $Id: icerhg_evp.F90 8378 2017-07-26 13:55:59Z clem $
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE ice_rhg_evp( stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
54      !!-------------------------------------------------------------------
55      !!                 ***  SUBROUTINE ice_rhg_evp  ***
56      !!                          EVP-C-grid
57      !!
58      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
59      !!  stress and sea-surface slope. Ice-ice interaction is described by
60      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
61      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
62      !!
63      !!  The points in the C-grid look like this, dear reader
64      !!
65      !!                              (ji,jj)
66      !!                                 |
67      !!                                 |
68      !!                      (ji-1,jj)  |  (ji,jj)
69      !!                             ---------   
70      !!                            |         |
71      !!                            | (ji,jj) |------(ji,jj)
72      !!                            |         |
73      !!                             ---------   
74      !!                     (ji-1,jj-1)     (ji,jj-1)
75      !!
76      !! ** Inputs  : - wind forcing (stress), oceanic currents
77      !!                ice total volume (vt_i) per unit area
78      !!                snow total volume (vt_s) per unit area
79      !!
80      !! ** Action  : - compute u_ice, v_ice : the components of the
81      !!                sea-ice velocity vector
82      !!              - compute delta_i, shear_i, divu_i, which are inputs
83      !!                of the ice thickness distribution
84      !!
85      !! ** Steps   : 1) Compute ice snow mass, ice strength
86      !!              2) Compute wind, oceanic stresses, mass terms and
87      !!                 coriolis terms of the momentum equation
88      !!              3) Solve the momentum equation (iterative procedure)
89      !!              4) Prevent high velocities if the ice is thin
90      !!              5) Recompute invariants of the strain rate tensor
91      !!                 which are inputs of the ITD, store stress
92      !!                 for the next time step
93      !!              6) Control prints of residual (convergence)
94      !!                 and charge ellipse.
95      !!                 The user should make sure that the parameters
96      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state
97      !!                 on the charge ellipse for plastic flow
98      !!                 e.g. in the Canadian Archipelago
99      !!
100      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013
101      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters)
102      !!              but this solution appears very unstable (see Kimmritz et al 2016)
103      !!
104      !! References : Hunke and Dukowicz, JPO97
105      !!              Bouillon et al., Ocean Modelling 2009
106      !!              Bouillon et al., Ocean Modelling 2013
107      !!-------------------------------------------------------------------
108      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   stress1_i, stress2_i, stress12_i
109      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   u_ice, v_ice, shear_i, divu_i, delta_i 
110      !!
111      INTEGER ::   ji, jj       ! dummy loop indices
112      INTEGER ::   jter         ! local integers
113      CHARACTER (len=50) ::   charout
114
115      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio
116      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
117      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
118      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
119      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
120      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
121      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars
122
123      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress
124      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
125      REAL(wp) ::   zintb, zintn                                             ! dummy argument
126      REAL(wp) ::   zfac_x, zfac_y
127     
128      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
129      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
130      !
131      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
132      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
133      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
134      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
135      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
136      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
137      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
138     
139      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
140      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
141      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
142      REAL(wp), DIMENSION(jpi,jpj) ::   zpice                           ! array used for the calculation of ice surface slope:
143                                                                             !   ocean surface (ssh_m) if ice is not embedded
144                                                                             !   ice top surface if ice is embedded   
145      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
146      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
147
148      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
149      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
150      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
151
152      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
153      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
154      !!-------------------------------------------------------------------
155
156#if defined key_agrif 
157      CALL agrif_interp_lim3( 'U', 0, nn_nevp )   ! First interpolation of coarse values
158      CALL agrif_interp_lim3( 'V', 0, nn_nevp )
159#endif
160      !
161      !------------------------------------------------------------------------------!
162      ! 0) mask at F points for the ice
163      !------------------------------------------------------------------------------!
164      ! ocean/land mask
165      DO jj = 1, jpjm1
166         DO ji = 1, jpim1      ! NO vector opt.
167            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
168         END DO
169      END DO
170      CALL lbc_lnk( zfmask, 'F', 1._wp )
171
172      ! Lateral boundary conditions on velocity (modify zfmask)
173      zwf(:,:) = zfmask(:,:)
174      DO jj = 2, jpjm1
175         DO ji = fs_2, fs_jpim1   ! vector opt.
176            IF( zfmask(ji,jj) == 0._wp ) THEN
177               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
178            ENDIF
179         END DO
180      END DO
181      DO jj = 2, jpjm1
182         IF( zfmask(1,jj) == 0._wp ) THEN
183            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
184         ENDIF
185         IF( zfmask(jpi,jj) == 0._wp ) THEN
186            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
187         ENDIF
188      END DO
189      DO ji = 2, jpim1
190         IF( zfmask(ji,1) == 0._wp ) THEN
191            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
192         ENDIF
193         IF( zfmask(ji,jpj) == 0._wp ) THEN
194            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
195         ENDIF
196      END DO
197      CALL lbc_lnk( zfmask, 'F', 1._wp )
198
199      !------------------------------------------------------------------------------!
200      ! 1) define some variables and initialize arrays
201      !------------------------------------------------------------------------------!
202      zrhoco = rau0 * rn_cio 
203
204      ! ecc2: square of yield ellipse eccenticrity
205      ecc2    = rn_ecc * rn_ecc
206      z1_ecc2 = 1._wp / ecc2
207
208      ! Time step for subcycling
209      zdtevp   = rdt_ice / REAL( nn_nevp )
210      z1_dtevp = 1._wp / zdtevp
211
212      ! alpha parameters (Bouillon 2009)
213      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
214      zalph2 = zalph1 * z1_ecc2
215
216      ! alpha and beta parameters (Bouillon 2013)
217      !!zalph1 = 40.
218      !!zalph2 = 40.
219      !!zbeta  = 3000.
220      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001)
221
222      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
223      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
224
225      ! Initialise stress tensor
226      zs1 (:,:) = stress1_i (:,:) 
227      zs2 (:,:) = stress2_i (:,:)
228      zs12(:,:) = stress12_i(:,:)
229
230      ! Ice strength
231      CALL ice_rdgrft_icestrength( nn_icestr )
232
233      ! scale factors
234      DO jj = 2, jpjm1
235         DO ji = fs_2, fs_jpim1
236            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
237            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
238         END DO
239      END DO
240           
241      !
242      !------------------------------------------------------------------------------!
243      ! 2) Wind / ocean stress, mass terms, coriolis terms
244      !------------------------------------------------------------------------------!
245
246      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
247         !                                           
248         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
249         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
250         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
251         !
252         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
253         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
254         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
255         !
256         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
257         !
258      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
259         zpice(:,:) = ssh_m(:,:)
260      ENDIF
261
262      DO jj = 2, jpjm1
263         DO ji = fs_2, fs_jpim1
264
265            ! ice fraction at U-V points
266            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
267            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
268
269            ! Ice/snow mass at U-V points
270            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
271            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
272            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
273            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
274            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
275
276            ! Ocean currents at U-V points
277            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
278               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
279           
280            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
281               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
282
283            ! Coriolis at T points (m*f)
284            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
285
286            ! m/dt
287            zmU_t(ji,jj)    = zmassU * z1_dtevp
288            zmV_t(ji,jj)    = zmassV * z1_dtevp
289
290            ! Drag ice-atm.
291            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
292            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
293
294            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
295            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
296            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
297
298            ! masks
299            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
300            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
301
302            ! switches
303            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
304            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
305
306         END DO
307      END DO
308      CALL lbc_lnk( zmf, 'T', 1. )
309      !
310      !------------------------------------------------------------------------------!
311      ! 3) Solution of the momentum equation, iterative procedure
312      !------------------------------------------------------------------------------!
313      !
314      !                                               !----------------------!
315      DO jter = 1 , nn_nevp                           !    loop over jter    !
316         !                                            !----------------------!       
317         IF(ln_ctl) THEN   ! Convergence test
318            DO jj = 1, jpjm1
319               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
320               zv_ice(:,jj) = v_ice(:,jj)
321            END DO
322         ENDIF
323
324         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
325         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
326            DO ji = 1, jpim1
327
328               ! shear at F points
329               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)   &
330                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
331                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
332
333            END DO
334         END DO
335         CALL lbc_lnk( zds, 'F', 1. )
336
337         DO jj = 2, jpjm1
338            DO ji = 2, jpim1 ! no vector loop
339
340               ! shear**2 at T points (doc eq. A16)
341               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
342                  &   + 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)  &
343                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
344             
345               ! divergence at T points
346               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
347                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
348                  &    ) * r1_e1e2t(ji,jj)
349               zdiv2 = zdiv * zdiv
350               
351               ! tension at T points
352               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)   &
353                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
354                  &   ) * r1_e1e2t(ji,jj)
355               zdt2 = zdt * zdt
356               
357               ! delta at T points
358               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
359
360               ! P/delta at T points
361               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
362               
363               ! stress at T points
364               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
365               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
366             
367            END DO
368         END DO
369         CALL lbc_lnk( zp_delt, 'T', 1. )
370
371         DO jj = 1, jpjm1
372            DO ji = 1, jpim1
373
374               ! P/delta at F points
375               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
376               
377               ! stress at F points
378               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
379
380            END DO
381         END DO
382         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
383 
384
385         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
386         DO jj = 2, jpjm1
387            DO ji = fs_2, fs_jpim1               
388               !                   !--- U points
389               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
390                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
391                  &                    ) * r1_e2u(ji,jj)                                                                      &
392                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
393                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
394                  &                  ) * r1_e1e2u(ji,jj)
395                  !
396                  !                !--- V points
397               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
398                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
399                  &                    ) * r1_e1v(ji,jj)                                                                      &
400                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
401                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
402                  &                  ) * r1_e1e2v(ji,jj)
403                  !
404                  !                !--- u_ice at V point
405               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
406                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
407                  !
408                  !                !--- v_ice at U point
409               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
410                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
411               !
412            END DO
413         END DO
414         !
415         ! --- Computation of ice velocity --- !
416         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
417         !  Bouillon et al. 2009 (eq 34-35) => stable
418         IF( MOD(jter,2) == 0 ) THEN ! even iterations
419            !
420            DO jj = 2, jpjm1
421               DO ji = fs_2, fs_jpim1
422                     !                 !--- tau_io/(v_oce - v_ice)
423                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
424                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
425                     !                 !--- Ocean-to-Ice stress
426                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
427                     !
428                     !                 !--- tau_bottom/v_ice
429                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
430                  zTauB = - tau_icebfr(ji,jj) / zvel
431                     !
432                     !                 !--- Coriolis at V-points (energy conserving formulation)
433                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
434                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
435                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
436                     !
437                     !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
438                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
439                     !
440                     !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
441                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
442                     !
443                     !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
444                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
445                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
446                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
447                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
448                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
449                     &           ) * zmaskV(ji,jj)
450                     !
451                  ! Bouillon 2013
452                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
453                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
454                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
455                  !
456               END DO
457            END DO
458            CALL lbc_lnk( v_ice, 'V', -1. )
459            !
460#if defined key_agrif
461!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
462            CALL agrif_interp_lim3( 'V' )
463#endif
464            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
465            !
466            DO jj = 2, jpjm1
467               DO ji = fs_2, fs_jpim1
468                               
469                  ! tau_io/(u_oce - u_ice)
470                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
471                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
472
473                  ! Ocean-to-Ice stress
474                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
475
476                  ! tau_bottom/u_ice
477                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
478                  zTauB = - tau_icebfr(ji,jj) / zvel
479
480                  ! Coriolis at U-points (energy conserving formulation)
481                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
482                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
483                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
484                 
485                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
486                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
487
488                  ! landfast switch => 0 = static friction ; 1 = sliding friction
489                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
490
491                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
492                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
493                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
494                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
495                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
496                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
497                     &           ) * zmaskU(ji,jj)
498                  ! Bouillon 2013
499                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
500                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
501                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
502               END DO
503            END DO
504            CALL lbc_lnk( u_ice, 'U', -1. )
505            !
506#if defined key_agrif
507!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
508            CALL agrif_interp_lim3( 'U' )
509#endif
510            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
511            !
512         ELSE ! odd iterations
513            !
514            DO jj = 2, jpjm1
515               DO ji = fs_2, fs_jpim1
516
517                  ! tau_io/(u_oce - u_ice)
518                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
519                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
520
521                  ! Ocean-to-Ice stress
522                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
523
524                  ! tau_bottom/u_ice
525                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
526                  zTauB = - tau_icebfr(ji,jj) / zvel
527
528                  ! Coriolis at U-points (energy conserving formulation)
529                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
530                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
531                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
532
533                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
534                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
535
536                  ! landfast switch => 0 = static friction ; 1 = sliding friction
537                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
538
539                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
540                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
541                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
542                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
543                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
544                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
545                     &           ) * zmaskU(ji,jj)
546                  ! Bouillon 2013
547                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
548                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
549                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
550               END DO
551            END DO
552            CALL lbc_lnk( u_ice, 'U', -1. )
553            !
554#if defined key_agrif
555!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
556            CALL agrif_interp_lim3( 'U' )
557#endif
558            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
559            !
560            DO jj = 2, jpjm1
561               DO ji = fs_2, fs_jpim1
562         
563                  ! tau_io/(v_oce - v_ice)
564                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
565                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
566
567                  ! Ocean-to-Ice stress
568                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
569
570                  ! tau_bottom/v_ice
571                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
572                  ztauB = - tau_icebfr(ji,jj) / zvel
573                 
574                  ! Coriolis at V-points (energy conserving formulation)
575                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
576                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
577                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
578
579                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
580                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
581
582                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
583                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
584                 
585                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
586                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
587                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
588                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
589                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
590                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
591                     &           ) * zmaskV(ji,jj)
592                  ! Bouillon 2013
593                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
594                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
595                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
596               END DO
597            END DO
598            CALL lbc_lnk( v_ice, 'V', -1. )
599            !
600#if defined key_agrif
601!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
602            CALL agrif_interp_lim3( 'V' )
603#endif
604            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
605            !
606         ENDIF
607         
608         IF(ln_ctl) THEN   ! Convergence test
609            DO jj = 2 , jpjm1
610               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
611            END DO
612            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
613            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
614         ENDIF
615         !
616         !                                                ! ==================== !
617      END DO                                              !  end loop over jter  !
618      !                                                   ! ==================== !
619      !
620      !------------------------------------------------------------------------------!
621      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
622      !------------------------------------------------------------------------------!
623      DO jj = 1, jpjm1
624         DO ji = 1, jpim1
625
626            ! shear at F points
627            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)   &
628               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
629               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
630
631         END DO
632      END DO           
633      CALL lbc_lnk( zds, 'F', 1. )
634     
635      DO jj = 2, jpjm1
636         DO ji = 2, jpim1 ! no vector loop
637           
638            ! tension**2 at T points
639            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)   &
640               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
641               &   ) * r1_e1e2t(ji,jj)
642            zdt2 = zdt * zdt
643           
644            ! shear**2 at T points (doc eq. A16)
645            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
646               &   + 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)  &
647               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
648           
649            ! shear at T points
650            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
651
652            ! divergence at T points
653            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
654               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
655               &            ) * r1_e1e2t(ji,jj)
656           
657            ! delta at T points
658            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
659            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
660            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
661
662         END DO
663      END DO
664      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
665     
666      ! --- Store the stress tensor for the next time step --- !
667      stress1_i (:,:) = zs1 (:,:)
668      stress2_i (:,:) = zs2 (:,:)
669      stress12_i(:,:) = zs12(:,:)
670      !
671
672      !------------------------------------------------------------------------------!
673      ! 5) SIMIP diagnostics
674      !------------------------------------------------------------------------------!
675
676!!gm  encapsulate with a flag (iom_use of the variable or better a flag defined one for all testing if one of the
677!!    diag is output.  NB the diag_...  are should only be ALLOCATED if the flag is true !
678
679      DO jj = 2, jpjm1
680         DO ji = 2, jpim1
681             rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
682
683             ! Stress tensor invariants (normal and shear stress N/m)
684             diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
685             diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
686
687             ! Stress terms of the momentum equation (N/m2)
688             diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term
689             diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
690
691             diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term
692             diag_corstry(ji,jj) = zCory(ji,jj) * rswitch
693
694             diag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term
695             diag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
696           
697             diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
698             diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
699
700             ! 2D ice mass, snow mass, area transport arrays (X, Y)
701             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
702             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
703
704             diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
705             diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
706
707             diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
708             diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
709
710             diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
711             diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
712
713         END DO
714      END DO
715
716      CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   &
717         &                  diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   &
718         &                  diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   & 
719         &                  diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    )
720
721      CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    )
722
723      CALL lbc_lnk_multi(   diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1.,   &
724         &                  diag_xatrp    , 'U', -1., diag_ymtrp_ice, 'V', -1.,   &
725         &                  diag_ymtrp_snw, 'V', -1., diag_yatrp    , 'V', -1.    )
726
727      !
728      !------------------------------------------------------------------------------!
729      ! 6) Control prints of residual and charge ellipse
730      !------------------------------------------------------------------------------!
731      !
732      ! print the residual for convergence
733      IF(ln_ctl) THEN
734         WRITE(charout,FMT="('ice_rhg_evp  : res =',D23.16, ' iter =',I4)") zresm, jter
735         CALL prt_ctl_info(charout)
736         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
737         !
738         ! print charge ellipse
739         ! This can be desactivated once the user is sure that the stress state
740      ! lie on the charge ellipse. See Bouillon et al. (2008) for more details
741         IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN
742            WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. '
743            CALL prt_ctl_info(charout)
744            DO jj = 2, jpjm1
745               DO ji = 2, jpim1
746                  IF (strength(ji,jj) > 1.0) THEN
747                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
748                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )
749                     WRITE(charout,FMT="('ice_rhg_evp  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
750                     CALL prt_ctl_info(charout)
751                  ENDIF
752               END DO
753            END DO
754            WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. '
755            CALL prt_ctl_info(charout)
756         ENDIF
757      ENDIF
758      !
759   END SUBROUTINE ice_rhg_evp
760
761#else
762   !!----------------------------------------------------------------------
763   !!   Default option         Empty module          NO LIM-3 sea-ice model
764   !!----------------------------------------------------------------------
765#endif
766
767   !!==============================================================================
768END MODULE icerhg_evp
Note: See TracBrowser for help on using the repository browser.