source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg_evp.F90 @ 11377

Last change on this file since 11377 was 11377, checked in by clem, 21 months ago

clean rheology, add a couple of outputs and remove landfast_home option

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