source: NEMO/branches/UKMO/dev_r10037_GPU/src/ICE/icedyn_rhg_evp.F90 @ 11467

Last change on this file since 11467 was 11467, checked in by andmirek, 22 months ago

Ticket #2197 allocate arrays at the beggining of the run

File size: 56.1 KB
Line 
1MODULE icedyn_rhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_evp  ***
4   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
5   !!======================================================================
6   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
7   !!            3.0  !  2008-03  (M. Vancoppenolle) adaptation to new model
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 case
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 + mEVP (Bouillon 2013)
13   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
14   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
15   !!----------------------------------------------------------------------
16#if defined key_si3
17   !!----------------------------------------------------------------------
18   !!   'key_si3'                                       SI3 sea-ice model
19   !!----------------------------------------------------------------------
20   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology
21   !!   rhg_evp_rst     : read/write EVP fields in ice restart
22   !!----------------------------------------------------------------------
23   USE phycst         ! Physical constant
24   USE dom_oce        ! Ocean domain
25   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
26   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
27   USE ice            ! sea-ice: ice variables
28   USE icedyn_rdgrft  ! sea-ice: ice strength
29   USE bdy_oce , ONLY : ln_bdy 
30   USE bdyice 
31#if defined key_agrif
32   USE agrif_ice_interp
33#endif
34   !
35   USE in_out_manager ! I/O manager
36   USE iom            ! I/O manager library
37   USE lib_mpp        ! MPP library
38   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
39   USE lbclnk         ! lateral boundary conditions (or mpp links)
40   USE prtctl         ! Print control
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   ice_dyn_rhg_evp   ! called by icedyn_rhg.F90
46   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
47
48   !! * Substitutions
49#  include "vectopt_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
58      !!-------------------------------------------------------------------
59      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
60      !!                             EVP-C-grid
61      !!
62      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
63      !!  stress and sea-surface slope. Ice-ice interaction is described by
64      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
65      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
66      !!
67      !!  The points in the C-grid look like this, dear reader
68      !!
69      !!                              (ji,jj)
70      !!                                 |
71      !!                                 |
72      !!                      (ji-1,jj)  |  (ji,jj)
73      !!                             ---------   
74      !!                            |         |
75      !!                            | (ji,jj) |------(ji,jj)
76      !!                            |         |
77      !!                             ---------   
78      !!                     (ji-1,jj-1)     (ji,jj-1)
79      !!
80      !! ** Inputs  : - wind forcing (stress), oceanic currents
81      !!                ice total volume (vt_i) per unit area
82      !!                snow total volume (vt_s) per unit area
83      !!
84      !! ** Action  : - compute u_ice, v_ice : the components of the
85      !!                sea-ice velocity vector
86      !!              - compute delta_i, shear_i, divu_i, which are inputs
87      !!                of the ice thickness distribution
88      !!
89      !! ** Steps   : 0) compute mask at F point
90      !!              1) Compute ice snow mass, ice strength
91      !!              2) Compute wind, oceanic stresses, mass terms and
92      !!                 coriolis terms of the momentum equation
93      !!              3) Solve the momentum equation (iterative procedure)
94      !!              4) Recompute delta, shear and divergence
95      !!                 (which are inputs of the ITD) & store stress
96      !!                 for the next time step
97      !!              5) Diagnostics including charge ellipse
98      !!
99      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
100      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
101      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
102      !!              (i.e. more stable and better convergence)
103      !!
104      !! References : Hunke and Dukowicz, JPO97
105      !!              Bouillon et al., Ocean Modelling 2009
106      !!              Bouillon et al., Ocean Modelling 2013
107      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
108      !!-------------------------------------------------------------------
109      USE scoce, ONLY : z1_e1t0 => scr2D1, z1_e2t0 => scr2D2, &                                       ! scale factors
110                        zp_delt => scr2D3, &                                                          ! P/delta at T points
111                        zbeta => scr2D4, &                                                            ! beta coef from Kimmritz 2017
112      !
113                        zdt_m => scr2D5, &                                                            ! (dt / ice-snow_mass) on T points
114                        zaU => scr2D6  , zaV => scr2D7, &                                             ! ice fraction on U/V points
115                        zmU_t => scr2D8, zmV_t => scr2D9, &                                           ! (ice-snow_mass / dt) on U/V points
116                        zmf => scr2D10, &                                                             ! coriolis parameter at T points
117                        zTauU_ia => scr2D11, ztauV_ia => scr2D12, &                                   ! ice-atm. stress at U-V points
118                        zspgU => scr2D13, zspgV => scr2D14, &                                         ! surface pressure gradient at U/V points
119                        v_oceU => scr2D15, u_oceV => scr2D16, v_iceU => scr2D17, u_iceV => scr2D18, & ! ocean/ice u/v component on V/U points                           
120                        zfU => scr2D19  , zfV => scr2D20, &                                           ! internal stresses
121      !
122                        zds => scr2D21, &                                                             ! shear
123                        zs1 => scr2D22, zs2 => scr2D23, zs12 => scr2D38, &                            ! stress tensor components
124                        zu_ice => scr2D24, zv_ice => scr2D25, zresr => scr2D26, &                     ! check convergence
125                        zpice  => scr2D27, &                                                          ! array used for the calculation of ice surface slope:
126      !                                                                                               !    ocean surface (ssh_m) if ice is not embedded
127      !                                                                                               !    ice top surface if ice is embedded   
128                        zCorx => scr2D28, zCory => scr2D29, &                                         ! Coriolis stress array
129                        ztaux_oi => scr2D30, ztauy_oi => scr2D31, &                                   ! Ocean-to-ice stress array
130      !
131                        zswitchU => scr2D32, zswitchV => scr2D33, &                                   ! dummy arrays
132                        zmaskU => scr2D34, zmaskV => scr2D35, &                                       ! mask for ice presence
133                        zfmask => scr2D36, zwf => scr2D37                                             ! mask at F points for the ice
134
135                       
136
137      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
138      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
139      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
140      !!
141      INTEGER ::   ji, jj       ! dummy loop indices
142      INTEGER ::   jter         ! local integers
143      !
144      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
145      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
146      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
147      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
148      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass
149      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
150      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars
151      !
152      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
153      REAL(wp) ::   zintb, zintn                                        ! dummy argument
154      REAL(wp) ::   zfac_x, zfac_y
155      REAL(wp) ::   zshear, zdum1, zdum2
156      !
157      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
158      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
159      !! --- diags
160      REAL(wp), DIMENSION(jpi,jpj) ::   zswi
161      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
162      !! --- SIMIP diags
163      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
164      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice
165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2)
166      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2)
167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2)
168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2)
169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2)
170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2)
171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress 
173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
175      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
178      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)     
179      !!-------------------------------------------------------------------
180
181      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
182      !
183!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
184      !------------------------------------------------------------------------------!
185      ! 0) mask at F points for the ice
186      !------------------------------------------------------------------------------!
187      ! ocean/land mask
188      DO jj = 1, jpjm1
189         DO ji = 1, jpim1      ! NO vector opt.
190            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
191         END DO
192      END DO
193      CALL lbc_lnk( zfmask, 'F', 1._wp )
194
195      ! Lateral boundary conditions on velocity (modify zfmask)
196      zwf(:,:) = zfmask(:,:)
197      DO jj = 2, jpjm1
198         DO ji = fs_2, fs_jpim1   ! vector opt.
199            IF( zfmask(ji,jj) == 0._wp ) THEN
200               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) ) )
201            ENDIF
202         END DO
203      END DO
204      DO jj = 2, jpjm1
205         IF( zfmask(1,jj) == 0._wp ) THEN
206            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
207         ENDIF
208         IF( zfmask(jpi,jj) == 0._wp ) THEN
209            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
210         ENDIF
211      END DO
212      DO ji = 2, jpim1
213         IF( zfmask(ji,1) == 0._wp ) THEN
214            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
215         ENDIF
216         IF( zfmask(ji,jpj) == 0._wp ) THEN
217            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
218         ENDIF
219      END DO
220      CALL lbc_lnk( zfmask, 'F', 1._wp )
221
222      !------------------------------------------------------------------------------!
223      ! 1) define some variables and initialize arrays
224      !------------------------------------------------------------------------------!
225      zrhoco = rau0 * rn_cio 
226
227      ! ecc2: square of yield ellipse eccenticrity
228      ecc2    = rn_ecc * rn_ecc
229      z1_ecc2 = 1._wp / ecc2
230
231      ! Time step for subcycling
232      zdtevp   = rdt_ice / REAL( nn_nevp )
233      z1_dtevp = 1._wp / zdtevp
234
235      ! alpha parameters (Bouillon 2009)
236      IF( .NOT. ln_aEVP ) THEN
237         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
238         zalph2 = zalph1 * z1_ecc2
239
240         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
241         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
242      ENDIF
243         
244      ! Initialise stress tensor
245      zs1 (:,:) = pstress1_i (:,:) 
246      zs2 (:,:) = pstress2_i (:,:)
247      zs12(:,:) = pstress12_i(:,:)
248
249      ! Ice strength
250      CALL ice_strength
251
252      ! scale factors
253      DO jj = 2, jpjm1
254         DO ji = fs_2, fs_jpim1
255            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
256            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
257         END DO
258      END DO
259           
260      !
261      !------------------------------------------------------------------------------!
262      ! 2) Wind / ocean stress, mass terms, coriolis terms
263      !------------------------------------------------------------------------------!
264
265      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
266         !                                           
267         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
268         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
269         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
270         !
271         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
272         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
273         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
274         !
275         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
276         !
277      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
278         zpice(:,:) = ssh_m(:,:)
279      ENDIF
280
281      DO jj = 2, jpjm1
282         DO ji = fs_2, fs_jpim1
283
284            ! ice fraction at U-V points
285            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)
286            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)
287
288            ! Ice/snow mass at U-V points
289            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
290            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
291            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
292            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
293            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
294
295            ! Ocean currents at U-V points
296            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
297               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
298           
299            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
300               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
301
302            ! Coriolis at T points (m*f)
303            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
304
305            ! dt/m at T points (for alpha and beta coefficients)
306            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
307           
308            ! m/dt
309            zmU_t(ji,jj)    = zmassU * z1_dtevp
310            zmV_t(ji,jj)    = zmassV * z1_dtevp
311           
312            ! Drag ice-atm.
313            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
314            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
315
316            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
317            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
318            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
319
320            ! masks
321            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
322            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
323
324            ! switches
325            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
326            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
327
328         END DO
329      END DO
330      CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. )
331      !
332      !------------------------------------------------------------------------------!
333      ! 3) Solution of the momentum equation, iterative procedure
334      !------------------------------------------------------------------------------!
335      !
336      !                                               !----------------------!
337      DO jter = 1 , nn_nevp                           !    loop over jter    !
338         !                                            !----------------------!       
339         IF(ln_ctl) THEN   ! Convergence test
340            DO jj = 1, jpjm1
341               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
342               zv_ice(:,jj) = v_ice(:,jj)
343            END DO
344         ENDIF
345
346         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
347         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
348            DO ji = 1, jpim1
349
350               ! shear at F points
351               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)   &
352                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
353                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
354
355            END DO
356         END DO
357         CALL lbc_lnk( zds, 'F', 1. )
358
359         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
360            DO ji = 2, jpi ! no vector loop
361
362               ! shear**2 at T points (doc eq. A16)
363               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
364                  &   + 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)  &
365                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
366             
367               ! divergence at T points
368               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
369                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
370                  &    ) * r1_e1e2t(ji,jj)
371               zdiv2 = zdiv * zdiv
372               
373               ! tension at T points
374               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)   &
375                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
376                  &   ) * r1_e1e2t(ji,jj)
377               zdt2 = zdt * zdt
378               
379               ! delta at T points
380               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
381
382               ! P/delta at T points
383               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
384
385               ! alpha & beta for aEVP
386               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
387               !   alpha = beta = sqrt(4*gamma)
388               IF( ln_aEVP ) THEN
389                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
390                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
391                  zalph2   = zalph1
392                  z1_alph2 = z1_alph1
393               ENDIF
394               
395               ! stress at T points
396               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
397               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
398             
399            END DO
400         END DO
401         CALL lbc_lnk( zp_delt, 'T', 1. )
402
403         DO jj = 1, jpjm1
404            DO ji = 1, jpim1
405
406               ! alpha & beta for aEVP
407               IF( ln_aEVP ) THEN
408                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
409                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
410                  zbeta(ji,jj) = zalph2
411               ENDIF
412               
413               ! P/delta at F points
414               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) )
415               
416               ! stress at F points
417               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
418
419            END DO
420         END DO
421
422         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
423         DO jj = 2, jpjm1
424            DO ji = fs_2, fs_jpim1               
425               !                   !--- U points
426               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
427                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
428                  &                    ) * r1_e2u(ji,jj)                                                                      &
429                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
430                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
431                  &                  ) * r1_e1e2u(ji,jj)
432               !
433               !                !--- V points
434               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
435                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
436                  &                    ) * r1_e1v(ji,jj)                                                                      &
437                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
438                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
439                  &                  ) * r1_e1e2v(ji,jj)
440               !
441               !                !--- u_ice at V point
442               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
443                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
444               !
445               !                !--- v_ice at U point
446               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
447                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
448            END DO
449         END DO
450         !
451         ! --- Computation of ice velocity --- !
452         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
453         !  Bouillon et al. 2009 (eq 34-35) => stable
454         IF( MOD(jter,2) == 0 ) THEN ! even iterations
455            !
456            DO jj = 2, jpjm1
457               DO ji = fs_2, fs_jpim1
458                  !                 !--- tau_io/(v_oce - v_ice)
459                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
460                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
461                  !                 !--- Ocean-to-Ice stress
462                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
463                  !
464                  !                 !--- tau_bottom/v_ice
465                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
466                  zTauB = - tau_icebfr(ji,jj) / zvel
467                  !
468                  !                 !--- Coriolis at V-points (energy conserving formulation)
469                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
470                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
471                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
472                  !
473                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
474                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
475                  !
476                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
477                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
478                  !
479                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
480                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
481                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
482                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
483                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
484                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
485                     &           ) * zmaskV(ji,jj)
486                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
487                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
488                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
489                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
490                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
491                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
492                     &            ) * zmaskV(ji,jj)
493                  ENDIF
494               END DO
495            END DO
496            CALL lbc_lnk( v_ice, 'V', -1. )
497            !
498#if defined key_agrif
499!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
500            CALL agrif_interp_ice( 'V' )
501#endif
502            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
503            !
504            DO jj = 2, jpjm1
505               DO ji = fs_2, fs_jpim1         
506                  !                 !--- tau_io/(u_oce - u_ice)
507                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
508                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
509                  !                 !--- Ocean-to-Ice stress
510                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
511                  !
512                  !                 !--- tau_bottom/u_ice
513                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
514                  zTauB = - tau_icebfr(ji,jj) / zvel
515                  !
516                  !                 !--- Coriolis at U-points (energy conserving formulation)
517                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
518                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
519                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
520                  !
521                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
522                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
523                  !
524                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
525                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
526                  !
527                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
528                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
529                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
530                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
531                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
532                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
533                     &            ) * zmaskU(ji,jj)
534                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
535                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
536                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
537                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
538                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
539                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
540                     &            ) * zmaskU(ji,jj)
541                  ENDIF
542               END DO
543            END DO
544            CALL lbc_lnk( u_ice, 'U', -1. )
545            !
546#if defined key_agrif
547!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
548            CALL agrif_interp_ice( 'U' )
549#endif
550            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
551            !
552         ELSE ! odd iterations
553            !
554            DO jj = 2, jpjm1
555               DO ji = fs_2, fs_jpim1
556                  !                 !--- tau_io/(u_oce - u_ice)
557                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
558                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
559                  !                 !--- Ocean-to-Ice stress
560                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
561                  !
562                  !                 !--- tau_bottom/u_ice
563                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
564                  zTauB = - tau_icebfr(ji,jj) / zvel
565                  !
566                  !                 !--- Coriolis at U-points (energy conserving formulation)
567                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
568                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
569                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
570                  !
571                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
572                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
573                  !
574                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
575                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
576                  !
577                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
578                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
579                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
580                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
581                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
582                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
583                     &            ) * zmaskU(ji,jj)
584                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
585                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
586                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
587                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
588                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
589                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
590                     &            ) * zmaskU(ji,jj)
591                  ENDIF
592               END DO
593            END DO
594            CALL lbc_lnk( u_ice, 'U', -1. )
595            !
596#if defined key_agrif
597!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
598            CALL agrif_interp_ice( 'U' )
599#endif
600            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
601            !
602            DO jj = 2, jpjm1
603               DO ji = fs_2, fs_jpim1
604                  !                 !--- tau_io/(v_oce - v_ice)
605                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
606                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
607                  !                 !--- Ocean-to-Ice stress
608                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
609                  !
610                  !                 !--- tau_bottom/v_ice
611                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
612                  ztauB = - tau_icebfr(ji,jj) / zvel
613                  !
614                  !                 !--- Coriolis at v-points (energy conserving formulation)
615                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
616                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
617                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
618                  !
619                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
620                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
621                  !
622                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
623                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
624                  !
625                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
626                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
627                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
628                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
629                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
630                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
631                     &           ) * zmaskV(ji,jj)
632                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
633                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
634                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
635                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
636                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
637                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
638                     &            ) * zmaskV(ji,jj)
639                  ENDIF
640               END DO
641            END DO
642            CALL lbc_lnk( v_ice, 'V', -1. )
643            !
644#if defined key_agrif
645!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
646            CALL agrif_interp_ice( 'V' )
647#endif
648            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
649            !
650         ENDIF
651         
652         IF(ln_ctl) THEN   ! Convergence test
653            DO jj = 2 , jpjm1
654               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
655            END DO
656            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
657            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
658         ENDIF
659         !
660         !                                                ! ==================== !
661      END DO                                              !  end loop over jter  !
662      !                                                   ! ==================== !
663      !
664      !------------------------------------------------------------------------------!
665      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
666      !------------------------------------------------------------------------------!
667      DO jj = 1, jpjm1
668         DO ji = 1, jpim1
669
670            ! shear at F points
671            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)   &
672               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
673               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
674
675         END DO
676      END DO           
677     
678      DO jj = 2, jpjm1
679         DO ji = 2, jpim1 ! no vector loop
680           
681            ! tension**2 at T points
682            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)   &
683               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
684               &   ) * r1_e1e2t(ji,jj)
685            zdt2 = zdt * zdt
686           
687            ! shear**2 at T points (doc eq. A16)
688            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
689               &   + 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)  &
690               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
691           
692            ! shear at T points
693            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
694
695            ! divergence at T points
696            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
697               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
698               &             ) * r1_e1e2t(ji,jj)
699           
700            ! delta at T points
701            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
702            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
703            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
704
705         END DO
706      END DO
707      CALL lbc_lnk_multi( pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
708     
709      ! --- Store the stress tensor for the next time step --- !
710      CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
711      pstress1_i (:,:) = zs1 (:,:)
712      pstress2_i (:,:) = zs2 (:,:)
713      pstress12_i(:,:) = zs12(:,:)
714      !
715
716      !------------------------------------------------------------------------------!
717      ! 5) diagnostics
718      !------------------------------------------------------------------------------!
719      DO jj = 1, jpj
720         DO ji = 1, jpi
721            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
722         END DO
723      END DO
724
725      ! --- divergence, shear and strength --- !
726      IF( iom_use('icediv') )   CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) )   ! divergence
727      IF( iom_use('iceshe') )   CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:) )   ! shear
728      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
729
730      ! --- charge ellipse --- !
731      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
732         !
733         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
734         !         
735         DO jj = 2, jpjm1
736            DO ji = 2, jpim1
737               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
738                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
739                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
740
741               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 
742
743               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
744
745!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
746!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
747!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
748!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
749               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
750               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
751               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
752            END DO
753         END DO
754         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
755         !
756         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
757         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
758         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
759         !
760         DEALLOCATE( zsig1 , zsig2 , zsig3 )
761      ENDIF
762     
763      ! --- SIMIP --- !
764      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
765         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
766         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
767         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
768
769         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
770            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
771            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
772            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
773         
774         DO jj = 2, jpjm1
775            DO ji = 2, jpim1
776               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
777               
778               ! Stress tensor invariants (normal and shear stress N/m)
779               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
780               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
781               
782               ! Stress terms of the momentum equation (N/m2)
783               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch     ! sea surface slope stress term
784               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
785               
786               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch     ! Coriolis stress term
787               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
788               
789               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch     ! internal stress term
790               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
791               
792               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
793               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
794               
795               ! 2D ice mass, snow mass, area transport arrays (X, Y)
796               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
797               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
798               
799               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
800               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
801               
802               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
803               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
804               
805               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
806               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
807               
808            END DO
809         END DO
810         
811         CALL lbc_lnk_multi( zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
812            &                zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
813            &                zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
814            &                zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
815                 
816         CALL lbc_lnk_multi( zdiag_utau_oi  , 'U', -1., zdiag_vtau_oi  , 'V', -1.,   &
817            &                zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
818            &                zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
819            &                zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
820         
821         IF( iom_use('normstr' ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
822         IF( iom_use('sheastr' ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
823         IF( iom_use('dssh_dx' ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
824         IF( iom_use('dssh_dy' ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
825         IF( iom_use('corstrx' ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
826         IF( iom_use('corstry' ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
827         IF( iom_use('intstrx' ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
828         IF( iom_use('intstry' ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
829         IF( iom_use('utau_oi' ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
830         IF( iom_use('vtau_oi' ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
831         IF( iom_use('xmtrpice') )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
832         IF( iom_use('ymtrpice') )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport
833         IF( iom_use('xmtrpsnw') )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
834         IF( iom_use('ymtrpsnw') )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
835         IF( iom_use('xatrp'   ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
836         IF( iom_use('yatrp'   ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
837
838         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
839            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
840            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
841            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
842
843      ENDIF
844      !
845   END SUBROUTINE ice_dyn_rhg_evp
846
847
848   SUBROUTINE rhg_evp_rst( cdrw, kt )
849      !!---------------------------------------------------------------------
850      !!                   ***  ROUTINE rhg_evp_rst  ***
851      !!                     
852      !! ** Purpose :   Read or write RHG file in restart file
853      !!
854      !! ** Method  :   use of IOM library
855      !!----------------------------------------------------------------------
856      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
857      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
858      !
859      INTEGER  ::   iter            ! local integer
860      INTEGER  ::   id1, id2, id3   ! local integers
861      !!----------------------------------------------------------------------
862      !
863      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
864         !                                   ! ---------------
865         IF( ln_rstart ) THEN                   !* Read the restart file
866            !
867            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
868            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
869            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
870            !
871            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
872               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
873               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
874               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
875            ELSE                                     ! start rheology from rest
876               IF(lwp) WRITE(numout,*)
877               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
878               stress1_i (:,:) = 0._wp
879               stress2_i (:,:) = 0._wp
880               stress12_i(:,:) = 0._wp
881            ENDIF
882         ELSE                                   !* Start from rest
883            IF(lwp) WRITE(numout,*)
884            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
885            stress1_i (:,:) = 0._wp
886            stress2_i (:,:) = 0._wp
887            stress12_i(:,:) = 0._wp
888         ENDIF
889         !
890      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
891         !                                   ! -------------------
892         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
893         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
894         !
895         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
896         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
897         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
898         !
899      ENDIF
900      !
901   END SUBROUTINE rhg_evp_rst
902
903#else
904   !!----------------------------------------------------------------------
905   !!   Default option         Empty module           NO SI3 sea-ice model
906   !!----------------------------------------------------------------------
907#endif
908
909   !!==============================================================================
910END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.