source: NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 @ 13286

Last change on this file since 13286 was 13286, checked in by smasson, 3 months ago

trunk: merge extra halos branch in trunk, see #2366

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