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

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

change calls in icestp.F90 for advection

File size: 41.4 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/LIM3 4.0 , UCL - NEMO Consortium (2011)
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
389               ! U points
390               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
391                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
392                  &                    ) * r1_e2u(ji,jj)                                                                      &
393                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
394                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
395                  &                  ) * r1_e1e2u(ji,jj)
396
397               ! V points
398               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
399                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
400                  &                    ) * r1_e1v(ji,jj)                                                                      &
401                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
402                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
403                  &                  ) * r1_e1e2v(ji,jj)
404
405               ! u_ice at V point
406               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
407                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
408               
409               ! v_ice at U point
410               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
411                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
412
413            END DO
414         END DO
415         !
416         ! --- Computation of ice velocity --- !
417         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
418         !  Bouillon et al. 2009 (eq 34-35) => stable
419         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
420           
421            DO jj = 2, jpjm1
422               DO ji = fs_2, fs_jpim1
423
424                  ! tau_io/(v_oce - v_ice)
425                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
426                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
427
428                  ! Ocean-to-Ice stress
429                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
430
431                  ! tau_bottom/v_ice
432                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
433                  zTauB = - tau_icebfr(ji,jj) / zvel
434
435                  ! Coriolis at V-points (energy conserving formulation)
436                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
437                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
438                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
439
440                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
441                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
442
443                  ! landfast switch => 0 = static friction ; 1 = sliding friction
444                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
445                 
446                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
447                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
448                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
449                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
450                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
451                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
452                     &           ) * zmaskV(ji,jj)
453                  ! Bouillon 2013
454                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
455                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
456                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
457
458               END DO
459            END DO
460            CALL lbc_lnk( v_ice, 'V', -1. )
461           
462#if defined key_agrif
463!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
464            CALL agrif_interp_lim3( 'V' )
465#endif
466            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
467
468            DO jj = 2, jpjm1
469               DO ji = fs_2, fs_jpim1
470                               
471                  ! tau_io/(u_oce - u_ice)
472                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
473                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
474
475                  ! Ocean-to-Ice stress
476                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
477
478                  ! tau_bottom/u_ice
479                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
480                  zTauB = - tau_icebfr(ji,jj) / zvel
481
482                  ! Coriolis at U-points (energy conserving formulation)
483                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
484                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
485                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
486                 
487                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
488                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
489
490                  ! landfast switch => 0 = static friction ; 1 = sliding friction
491                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
492
493                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
494                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
495                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
496                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
497                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
498                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
499                     &           ) * zmaskU(ji,jj)
500                  ! Bouillon 2013
501                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
502                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
503                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
504               END DO
505            END DO
506            CALL lbc_lnk( u_ice, 'U', -1. )
507           
508#if defined key_agrif
509!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
510            CALL agrif_interp_lim3( 'U' )
511#endif
512            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
513
514         ELSE ! odd iterations
515
516            DO jj = 2, jpjm1
517               DO ji = fs_2, fs_jpim1
518
519                  ! tau_io/(u_oce - u_ice)
520                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
521                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
522
523                  ! Ocean-to-Ice stress
524                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
525
526                  ! tau_bottom/u_ice
527                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
528                  zTauB = - tau_icebfr(ji,jj) / zvel
529
530                  ! Coriolis at U-points (energy conserving formulation)
531                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
532                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
533                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
534
535                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
536                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
537
538                  ! landfast switch => 0 = static friction ; 1 = sliding friction
539                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
540
541                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
542                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
543                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
544                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
545                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
546                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
547                     &           ) * zmaskU(ji,jj)
548                  ! Bouillon 2013
549                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
550                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
551                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
552               END DO
553            END DO
554            CALL lbc_lnk( u_ice, 'U', -1. )
555           
556#if defined key_agrif
557!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
558            CALL agrif_interp_lim3( 'U' )
559#endif
560            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
561
562            DO jj = 2, jpjm1
563               DO ji = fs_2, fs_jpim1
564         
565                  ! tau_io/(v_oce - v_ice)
566                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
567                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
568
569                  ! Ocean-to-Ice stress
570                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
571
572                  ! tau_bottom/v_ice
573                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
574                  ztauB = - tau_icebfr(ji,jj) / zvel
575                 
576                  ! Coriolis at V-points (energy conserving formulation)
577                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
578                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
579                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
580
581                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
582                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
583
584                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
585                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
586                 
587                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
588                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
589                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
590                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
591                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
592                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
593                     &           ) * zmaskV(ji,jj)
594                  ! Bouillon 2013
595                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
596                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
597                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
598               END DO
599            END DO
600            CALL lbc_lnk( v_ice, 'V', -1. )
601           
602#if defined key_agrif
603!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
604            CALL agrif_interp_lim3( 'V' )
605#endif
606            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
607
608         ENDIF
609         
610         IF(ln_ctl) THEN   ! Convergence test
611            DO jj = 2 , jpjm1
612               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
613            END DO
614            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
615            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
616         ENDIF
617         !
618         !                                                ! ==================== !
619      END DO                                              !  end loop over jter  !
620      !                                                   ! ==================== !
621      !
622      !------------------------------------------------------------------------------!
623      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
624      !------------------------------------------------------------------------------!
625      DO jj = 1, jpjm1
626         DO ji = 1, jpim1
627
628            ! shear at F points
629            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)   &
630               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
631               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
632
633         END DO
634      END DO           
635      CALL lbc_lnk( zds, 'F', 1. )
636     
637      DO jj = 2, jpjm1
638         DO ji = 2, jpim1 ! no vector loop
639           
640            ! tension**2 at T points
641            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)   &
642               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
643               &   ) * r1_e1e2t(ji,jj)
644            zdt2 = zdt * zdt
645           
646            ! shear**2 at T points (doc eq. A16)
647            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
648               &   + 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)  &
649               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
650           
651            ! shear at T points
652            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
653
654            ! divergence at T points
655            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
656               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
657               &            ) * r1_e1e2t(ji,jj)
658           
659            ! delta at T points
660            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
661            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
662            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
663
664         END DO
665      END DO
666      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
667     
668      ! --- Store the stress tensor for the next time step --- !
669      stress1_i (:,:) = zs1 (:,:)
670      stress2_i (:,:) = zs2 (:,:)
671      stress12_i(:,:) = zs12(:,:)
672      !
673
674      !------------------------------------------------------------------------------!
675      ! 5) SIMIP diagnostics
676      !------------------------------------------------------------------------------!
677                           
678      DO jj = 2, jpjm1
679         DO ji = 2, jpim1
680             rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
681
682             ! Stress tensor invariants (normal and shear stress N/m)
683             diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
684             diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
685
686             ! Stress terms of the momentum equation (N/m2)
687             diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term
688             diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
689
690             diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term
691             diag_corstry(ji,jj) = zCory(ji,jj) * rswitch
692
693             diag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term
694             diag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
695           
696             diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
697             diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
698
699             ! 2D ice mass, snow mass, area transport arrays (X, Y)
700             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
701             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
702
703             diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
704             diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
705
706             diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
707             diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
708
709             diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
710             diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
711
712         END DO
713      END DO
714
715      CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   &
716                 &          diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   &
717                 &          diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   & 
718                 &          diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    )
719
720      CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    )
721
722      CALL lbc_lnk_multi(   diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1., &
723                 &          diag_xatrp    , 'U', -1., diag_ymtrp_ice, 'V', -1., &
724                 &          diag_ymtrp_snw, 'V', -1., diag_yatrp    , 'V', -1.  )
725
726      !
727      !------------------------------------------------------------------------------!
728      ! 6) Control prints of residual and charge ellipse
729      !------------------------------------------------------------------------------!
730      !
731      ! print the residual for convergence
732      IF(ln_ctl) THEN
733         WRITE(charout,FMT="('ice_rhg_evp  : res =',D23.16, ' iter =',I4)") zresm, jter
734         CALL prt_ctl_info(charout)
735         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
736      ENDIF
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. 08 for more details
741      IF(ln_ctl) THEN
742         IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN
743            WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. '
744            CALL prt_ctl_info(charout)
745            DO jj = 2, jpjm1
746               DO ji = 2, jpim1
747                  IF (strength(ji,jj) > 1.0) THEN
748                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
749                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )
750                     WRITE(charout,FMT="('ice_rhg_evp  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
751                     CALL prt_ctl_info(charout)
752                  ENDIF
753               END DO
754            END DO
755            WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. '
756            CALL prt_ctl_info(charout)
757         ENDIF
758      ENDIF
759      !
760   END SUBROUTINE ice_rhg_evp
761
762#endif
763
764   !!==============================================================================
765END MODULE icerhg_evp
Note: See TracBrowser for help on using the repository browser.