source: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90 @ 11877

Last change on this file since 11877 was 11877, checked in by stefryn, 11 months ago

add initialisation and restart of sturcture tensor

File size: 56.3 KB
Line 
1MODULE icedyn_rhg_eap
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_eap  ***
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_eap : computes ice velocities from EVP rheology
21   !!   rhg_eap_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_eap   ! called by icedyn_rhg.F90
47   PUBLIC   rhg_eap_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: icedyn_rhg_eap.F90 11536 2019-09-11 13:54:18Z smasson $
54   !! Software governed by the CeCILL license (see ./LICENSE)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, aniso_11, aniso_12 )
59      !!-------------------------------------------------------------------
60      !!                 ***  SUBROUTINE ice_dyn_rhg_eap  ***
61      !!                             EAP-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-anisotropic-plastic (EAP) law including shear
66      !!  strength and a bulk rheology . 
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      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   aniso_11, aniso_12   ! structure tensor components
114      !!
115      INTEGER ::   ji, jj       ! dummy loop indices
116      INTEGER ::   jter         ! local integers
117      !
118      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
119      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
120      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
121      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
122      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
123      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
124      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
125      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
126      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
127      !
128      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
129      REAL(wp) ::   zintb, zintn                                        ! dummy argument
130      REAL(wp) ::   zfac_x, zfac_y
131      REAL(wp) ::   zshear, zdum1, zdum2
132      !
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_eap: EAP 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_eap', 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_eap', 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      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
247      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile
248      ELSE                         ;   zkt = 0._wp
249      ENDIF
250      !
251      !------------------------------------------------------------------------------!
252      ! 2) Wind / ocean stress, mass terms, coriolis terms
253      !------------------------------------------------------------------------------!
254      ! sea surface height
255      !    embedded sea ice: compute representative ice top surface
256      !    non-embedded sea ice: use ocean surface for slope calculation
257      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
258
259      DO jj = 2, jpjm1
260         DO ji = fs_2, fs_jpim1
261
262            ! ice fraction at U-V points
263            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)
264            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)
265
266            ! Ice/snow mass at U-V points
267            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
268            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
269            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
270            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
271            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
272
273            ! Ocean currents at U-V points
274            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)
275            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)
276
277            ! Coriolis at T points (m*f)
278            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
279
280            ! dt/m at T points (for alpha and beta coefficients)
281            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
282           
283            ! m/dt
284            zmU_t(ji,jj)    = zmassU * z1_dtevp
285            zmV_t(ji,jj)    = zmassV * z1_dtevp
286           
287            ! Drag ice-atm.
288            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
289            ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
290
291            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
292            zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
293            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
294
295            ! masks
296            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
297            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
298
299            ! switches
300            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
301            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
302            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
303            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
304
305         END DO
306      END DO
307      CALL lbc_lnk_multi( 'icedyn_rhg_eap', zmf, 'T', 1., zdt_m, 'T', 1. )
308      !
309      !                                  !== Landfast ice parameterization ==!
310      !
311      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
312         DO jj = 2, jpjm1
313            DO ji = fs_2, fs_jpim1
314               ! ice thickness at U-V points
315               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)
316               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)
317               ! ice-bottom stress at U points
318               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
319               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
320               ! ice-bottom stress at V points
321               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
322               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
323               ! ice_bottom stress at T points
324               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
325               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
326            END DO
327         END DO
328         CALL lbc_lnk( 'icedyn_rhg_eap', tau_icebfr(:,:), 'T', 1. )
329         !
330      ELSE                               !-- no landfast
331         DO jj = 2, jpjm1
332            DO ji = fs_2, fs_jpim1
333               ztaux_base(ji,jj) = 0._wp
334               ztauy_base(ji,jj) = 0._wp
335            END DO
336         END DO
337      ENDIF
338
339      !------------------------------------------------------------------------------!
340      ! 3) Solution of the momentum equation, iterative procedure
341      !------------------------------------------------------------------------------!
342      !
343      !                                               ! ==================== !
344      DO jter = 1 , nn_nevp                           !    loop over jter    !
345         !                                            ! ==================== !       
346         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
347         !
348!!$         IF(ln_ctl) THEN   ! Convergence test
349!!$            DO jj = 1, jpjm1
350!!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
351!!$               zv_ice(:,jj) = v_ice(:,jj)
352!!$            END DO
353!!$         ENDIF
354
355         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
356         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
357            DO ji = 1, jpim1
358
359               ! shear at F points
360               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)   &
361                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
362                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
363
364            END DO
365         END DO
366         CALL lbc_lnk( 'icedyn_rhg_eap', zds, 'F', 1. )
367
368         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
369            DO ji = 2, jpi ! no vector loop
370
371               ! shear**2 at T points (doc eq. A16)
372               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
373                  &   + 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)  &
374                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
375             
376               ! divergence at T points
377               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
378                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
379                  &    ) * r1_e1e2t(ji,jj)
380               zdiv2 = zdiv * zdiv
381               
382               ! tension at T points
383               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)   &
384                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
385                  &   ) * r1_e1e2t(ji,jj)
386               zdt2 = zdt * zdt
387               
388               ! delta at T points
389               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
390
391               ! P/delta at T points
392               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
393
394               ! alpha & beta for aEVP
395               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
396               !   alpha = beta = sqrt(4*gamma)
397               IF( ln_aEVP ) THEN
398                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
399                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
400                  zalph2   = zalph1
401                  z1_alph2 = z1_alph1
402               ENDIF
403               
404               ! stress at T points (zkt/=0 if landfast)
405               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
406               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
407             
408            END DO
409         END DO
410         CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. )
411
412         DO jj = 1, jpjm1
413            DO ji = 1, jpim1
414
415               ! alpha & beta for aEVP
416               IF( ln_aEVP ) THEN
417                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
418                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
419                  zbeta(ji,jj) = zalph2
420               ENDIF
421               
422               ! P/delta at F points
423               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) )
424               
425               ! stress at F points (zkt/=0 if landfast)
426               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
427
428            END DO
429         END DO
430
431         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
432         DO jj = 2, jpjm1
433            DO ji = fs_2, fs_jpim1               
434               !                   !--- U points
435               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
436                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
437                  &                    ) * r1_e2u(ji,jj)                                                                      &
438                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
439                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
440                  &                  ) * r1_e1e2u(ji,jj)
441               !
442               !                !--- V points
443               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
444                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
445                  &                    ) * r1_e1v(ji,jj)                                                                      &
446                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
447                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
448                  &                  ) * r1_e1e2v(ji,jj)
449               !
450               !                !--- ice currents at U-V point
451               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)
452               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)
453               !
454            END DO
455         END DO
456         !
457         ! --- Computation of ice velocity --- !
458         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
459         !  Bouillon et al. 2009 (eq 34-35) => stable
460         IF( MOD(jter,2) == 0 ) THEN ! even iterations
461            !
462            DO jj = 2, jpjm1
463               DO ji = fs_2, fs_jpim1
464                  !                 !--- tau_io/(v_oce - v_ice)
465                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
466                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
467                  !                 !--- Ocean-to-Ice stress
468                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
469                  !
470                  !                 !--- tau_bottom/v_ice
471                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
472                  zTauB = ztauy_base(ji,jj) / zvel
473                  !                 !--- OceanBottom-to-Ice stress
474                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
475                  !
476                  !                 !--- Coriolis at V-points (energy conserving formulation)
477                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
478                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
479                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
480                  !
481                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
482                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
483                  !
484                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
485                  !                                         1 = sliding friction : TauB < RHS
486                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
487                  !
488                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
489                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
490                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
491                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
492                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
493                        &             ) * 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
494                        &           )   * zmsk00y(ji,jj)
495                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
496                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
497                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
498                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
499                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
500                        &              ) * 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
501                        &            )   * zmsk00y(ji,jj)
502                  ENDIF
503               END DO
504            END DO
505            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1. )
506            !
507#if defined key_agrif
508!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
509            CALL agrif_interp_ice( 'V' )
510#endif
511            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
512            !
513            DO jj = 2, jpjm1
514               DO ji = fs_2, fs_jpim1         
515                  !                 !--- tau_io/(u_oce - u_ice)
516                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
517                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
518                  !                 !--- Ocean-to-Ice stress
519                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
520                  !
521                  !                 !--- tau_bottom/u_ice
522                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
523                  zTauB = ztaux_base(ji,jj) / zvel
524                  !                 !--- OceanBottom-to-Ice stress
525                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
526                  !
527                  !                 !--- Coriolis at U-points (energy conserving formulation)
528                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
529                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
530                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
531                  !
532                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
533                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
534                  !
535                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
536                  !                                         1 = sliding friction : TauB < RHS
537                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
538                  !
539                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
540                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
541                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
542                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
543                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
544                        &             ) * 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
545                        &           )   * zmsk00x(ji,jj)
546                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
547                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
548                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
549                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
550                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
551                        &              ) * 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
552                        &            )   * zmsk00x(ji,jj)
553                  ENDIF
554               END DO
555            END DO
556            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1. )
557            !
558#if defined key_agrif
559!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
560            CALL agrif_interp_ice( 'U' )
561#endif
562            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
563            !
564         ELSE ! odd iterations
565            !
566            DO jj = 2, jpjm1
567               DO ji = fs_2, fs_jpim1
568                  !                 !--- tau_io/(u_oce - u_ice)
569                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
570                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
571                  !                 !--- Ocean-to-Ice stress
572                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
573                  !
574                  !                 !--- tau_bottom/u_ice
575                  zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
576                  zTauB = ztaux_base(ji,jj) / zvel
577                  !                 !--- OceanBottom-to-Ice stress
578                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
579                  !
580                  !                 !--- Coriolis at U-points (energy conserving formulation)
581                  zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
582                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
583                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
584                  !
585                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
586                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
587                  !
588                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
589                  !                                         1 = sliding friction : TauB < RHS
590                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
591                  !
592                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
593                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
594                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
595                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
596                        &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
597                        &             ) * 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
598                        &           )   * zmsk00x(ji,jj)
599                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
600                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
601                        &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
602                        &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
603                        &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
604                        &              ) * 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
605                        &            )   * zmsk00x(ji,jj)
606                  ENDIF
607               END DO
608            END DO
609            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1. )
610            !
611#if defined key_agrif
612!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
613            CALL agrif_interp_ice( 'U' )
614#endif
615            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
616            !
617            DO jj = 2, jpjm1
618               DO ji = fs_2, fs_jpim1
619                  !                 !--- tau_io/(v_oce - v_ice)
620                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
621                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
622                  !                 !--- Ocean-to-Ice stress
623                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
624                  !
625                  !                 !--- tau_bottom/v_ice
626                  zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
627                  zTauB = ztauy_base(ji,jj) / zvel
628                  !                 !--- OceanBottom-to-Ice stress
629                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
630                  !
631                  !                 !--- Coriolis at v-points (energy conserving formulation)
632                  zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
633                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
634                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
635                  !
636                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
637                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
638                  !
639                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
640                  !                                         1 = sliding friction : TauB < RHS
641                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
642                  !
643                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
644                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
645                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
646                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
647                        &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
648                        &             ) * 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
649                        &           )   * zmsk00y(ji,jj)
650                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
651                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
652                        &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
653                        &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast
654                        &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0
655                        &              ) * 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
656                        &            )   * zmsk00y(ji,jj)
657                  ENDIF
658               END DO
659            END DO
660            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1. )
661            !
662#if defined key_agrif
663!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
664            CALL agrif_interp_ice( 'V' )
665#endif
666            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
667            !
668         ENDIF
669
670!!$         IF(ln_ctl) THEN   ! Convergence test
671!!$            DO jj = 2 , jpjm1
672!!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
673!!$            END DO
674!!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
675!!$            CALL mpp_max( 'icedyn_rhg_eap', zresm )   ! max over the global domain
676!!$         ENDIF
677         !
678         !                                                ! ==================== !
679      END DO                                              !  end loop over jter  !
680      !                                                   ! ==================== !
681      !
682      !------------------------------------------------------------------------------!
683      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
684      !------------------------------------------------------------------------------!
685      DO jj = 1, jpjm1
686         DO ji = 1, jpim1
687
688            ! shear at F points
689            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)   &
690               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
691               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
692
693         END DO
694      END DO           
695     
696      DO jj = 2, jpjm1
697         DO ji = 2, jpim1 ! no vector loop
698           
699            ! tension**2 at T points
700            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)   &
701               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
702               &   ) * r1_e1e2t(ji,jj)
703            zdt2 = zdt * zdt
704           
705            ! shear**2 at T points (doc eq. A16)
706            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
707               &   + 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)  &
708               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
709           
710            ! shear at T points
711            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
712
713            ! divergence at T points
714            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
715               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
716               &             ) * r1_e1e2t(ji,jj)
717           
718            ! delta at T points
719            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
720            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
721            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
722
723         END DO
724      END DO
725      CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
726     
727      ! --- Store the stress tensor for the next time step --- !
728      CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
729      pstress1_i (:,:) = zs1 (:,:)
730      pstress2_i (:,:) = zs2 (:,:)
731      pstress12_i(:,:) = zs12(:,:)
732      !
733
734      !------------------------------------------------------------------------------!
735      ! 5) diagnostics
736      !------------------------------------------------------------------------------!
737      DO jj = 1, jpj
738         DO ji = 1, jpi
739            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
740         END DO
741      END DO
742
743      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
744      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
745         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
746         !
747         CALL lbc_lnk_multi( 'icedyn_rhg_eap', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
748            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
749         !
750         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
751         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
752         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
753         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
754         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
755         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
756      ENDIF
757       
758      ! --- divergence, shear and strength --- !
759      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
760      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
761      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
762
763      ! --- stress tensor --- !
764      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN
765         !
766         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
767         !         
768         DO jj = 2, jpjm1
769            DO ji = 2, jpim1
770               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
771                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
772                  &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )
773
774               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 
775
776               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )
777
778!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
779!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
780!!               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
781!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
782               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
783               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
784               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
785            END DO
786         END DO
787         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
788         !
789         CALL iom_put( 'isig1' , zsig1 )
790         CALL iom_put( 'isig2' , zsig2 )
791         CALL iom_put( 'isig3' , zsig3 )
792         !
793         ! Stress tensor invariants (normal and shear stress N/m)
794         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress
795         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress
796
797         DEALLOCATE( zsig1 , zsig2 , zsig3 )
798      ENDIF
799     
800      ! --- SIMIP --- !
801      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
802         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
803         !
804         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zspgU, 'U', -1., zspgV, 'V', -1., &
805            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )
806
807         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
808         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
809         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
810         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
811         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
812         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
813      ENDIF
814
815      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
816         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
817         !
818         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
819            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
820         !
821         DO jj = 2, jpjm1
822            DO ji = 2, jpim1
823               ! 2D ice mass, snow mass, area transport arrays (X, Y)
824               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
825               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
826
827               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
828               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
829
830               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
831               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
832
833               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
834               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
835
836            END DO
837         END DO
838
839         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
840            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
841            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
842
843         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
844         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport
845         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
846         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
847         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
848         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
849
850         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
851            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
852
853      ENDIF
854      !
855   END SUBROUTINE ice_dyn_rhg_eap
856
857
858   SUBROUTINE rhg_eap_rst( cdrw, kt )
859      !!---------------------------------------------------------------------
860      !!                   ***  ROUTINE rhg_eap_rst  ***
861      !!                     
862      !! ** Purpose :   Read or write RHG file in restart file
863      !!
864      !! ** Method  :   use of IOM library
865      !!----------------------------------------------------------------------
866      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
867      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
868      !
869      INTEGER  ::   iter            ! local integer
870      INTEGER  ::   id1, id2, id3, id4, id5   ! local integers
871      !!----------------------------------------------------------------------
872      !
873      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
874         !                                   ! ---------------
875         IF( ln_rstart ) THEN                   !* Read the restart file
876            !
877            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
878            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
879            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
880            id4 = iom_varid( numrir, 'aniso_11' , ldstop = .FALSE. )
881            id5 = iom_varid( numrir, 'aniso_12' , ldstop = .FALSE. )
882            !
883            IF( MIN( id1, id2, id3, id4, id5 ) > 0 ) THEN      ! fields exist
884               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
885               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
886               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
887               CALL iom_get( numrir, jpdom_autoglo, 'aniso_11' , aniso_11 )
888               CALL iom_get( numrir, jpdom_autoglo, 'aniso_12' , aniso_12 )
889            ELSE                                     ! start rheology from rest
890               IF(lwp) WRITE(numout,*)
891               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
892               stress1_i (:,:) = 0._wp
893               stress2_i (:,:) = 0._wp
894               stress12_i(:,:) = 0._wp
895               aniso_11 (:,:) = 0.5_wp
896               aniso_12 (:,:) = 0._wp
897            ENDIF
898         ELSE                                   !* Start from rest
899            IF(lwp) WRITE(numout,*)
900            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
901            stress1_i (:,:) = 0._wp
902            stress2_i (:,:) = 0._wp
903            stress12_i(:,:) = 0._wp
904            aniso_11 (:,:) = 0.5_wp
905            aniso_12 (:,:) = 0._wp
906         ENDIF
907         !
908      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
909         !                                   ! -------------------
910         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
911         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
912         !
913         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
914         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
915         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
916         CALL iom_rstput( iter, nitrst, numriw, 'aniso_11' , aniso_11 )
917         CALL iom_rstput( iter, nitrst, numriw, 'aniso_12' , aniso_12 )
918         !
919      ENDIF
920      !
921   END SUBROUTINE rhg_eap_rst
922
923#else
924   !!----------------------------------------------------------------------
925   !!   Default option         Empty module           NO SI3 sea-ice model
926   !!----------------------------------------------------------------------
927#endif
928
929   !!==============================================================================
930END MODULE icedyn_rhg_eap
Note: See TracBrowser for help on using the repository browser.