New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
icedyn_rhg_eap.F90 in NEMO/trunk/src/ICE – NEMO

source: NEMO/trunk/src/ICE/icedyn_rhg_eap.F90 @ 15062

Last change on this file since 15062 was 15062, checked in by jchanut, 12 months ago

Suppress time varying scale factors and depths declarations with key_qco and key_linssh. Remove spaces that preclude from correct replacement of some scale factor arrays during preprocessing stage (at least with Apple clang version 11.0.3, this is problem).

File size: 93.2 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   !!                 !  2019     (S. Rynders, Y. Aksenov, C. Rousset)  change into eap rheology from
16   !!                                           CICE code (Tsamados, Heorton)
17   !!----------------------------------------------------------------------
18#if defined key_si3
19   !!----------------------------------------------------------------------
20   !!   'key_si3'                                       SI3 sea-ice model
21   !!----------------------------------------------------------------------
22   !!   ice_dyn_rhg_eap : computes ice velocities from EVP rheology
23   !!   rhg_eap_rst     : read/write EVP fields in ice restart
24   !!----------------------------------------------------------------------
25   USE phycst         ! Physical constant
26   USE dom_oce        ! Ocean domain
27   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
28   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
29   USE ice            ! sea-ice: ice variables
30   USE icevar         ! ice_var_sshdyn
31   USE icedyn_rdgrft  ! sea-ice: ice strength
32   USE bdy_oce , ONLY : ln_bdy
33   USE bdyice
34#if defined key_agrif
35   USE agrif_ice_interp
36#endif
37   !
38   USE in_out_manager ! I/O manager
39   USE iom            ! I/O manager library
40   USE lib_mpp        ! MPP library
41   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
42   USE lbclnk         ! lateral boundary conditions (or mpp links)
43   USE prtctl         ! Print control
44
45   USE netcdf         ! NetCDF library for convergence test
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC   ice_dyn_rhg_eap   ! called by icedyn_rhg.F90
50   PUBLIC   rhg_eap_rst       ! called by icedyn_rhg.F90
51
52   REAL(wp), PARAMETER ::   pphi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg)
53
54   ! look-up table for calculating structure tensor
55   INTEGER, PARAMETER ::   nx_yield = 41
56   INTEGER, PARAMETER ::   ny_yield = 41
57   INTEGER, PARAMETER ::   na_yield = 21
58
59   REAL(wp), DIMENSION(nx_yield, ny_yield, na_yield) ::   s11r, s12r, s22r, s11s, s12s, s22s
60   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   fimask   ! mask at F points for the ice
61
62   !! for convergence tests
63   INTEGER ::   ncvgid   ! netcdf file id
64   INTEGER ::   nvarid   ! netcdf variable id
65
66   !! * Substitutions
67#  include "do_loop_substitute.h90"
68#  include "domzgr_substitute.h90"
69   !!----------------------------------------------------------------------
70   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
71   !! $Id: icedyn_rhg_eap.F90 11536 2019-09-11 13:54:18Z smasson $
72   !! Software governed by the CeCILL license (see ./LICENSE)
73   !!----------------------------------------------------------------------
74CONTAINS
75
76   SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12, prdg_conv )
77      !!-------------------------------------------------------------------
78      !!                 ***  SUBROUTINE ice_dyn_rhg_eap  ***
79      !!                             EAP-C-grid
80      !!
81      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
82      !!  stress and sea-surface slope. Ice-ice interaction is described by
83      !!  a non-linear elasto-anisotropic-plastic (EAP) law including shear
84      !!  strength and a bulk rheology .
85      !!
86      !!  The points in the C-grid look like this, dear reader
87      !!
88      !!                              (ji,jj)
89      !!                                 |
90      !!                                 |
91      !!                      (ji-1,jj)  |  (ji,jj)
92      !!                             ---------
93      !!                            |         |
94      !!                            | (ji,jj) |------(ji,jj)
95      !!                            |         |
96      !!                             ---------
97      !!                     (ji-1,jj-1)     (ji,jj-1)
98      !!
99      !! ** Inputs  : - wind forcing (stress), oceanic currents
100      !!                ice total volume (vt_i) per unit area
101      !!                snow total volume (vt_s) per unit area
102      !!
103      !! ** Action  : - compute u_ice, v_ice : the components of the
104      !!                sea-ice velocity vector
105      !!              - compute delta_i, shear_i, divu_i, which are inputs
106      !!                of the ice thickness distribution
107      !!
108      !! ** Steps   : 0) compute mask at F point
109      !!              1) Compute ice snow mass, ice strength
110      !!              2) Compute wind, oceanic stresses, mass terms and
111      !!                 coriolis terms of the momentum equation
112      !!              3) Solve the momentum equation (iterative procedure)
113      !!              4) Recompute delta, shear and divergence
114      !!                 (which are inputs of the ITD) & store stress
115      !!                 for the next time step
116      !!              5) Diagnostics including charge ellipse
117      !!
118      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
119      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
120      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
121      !!              (i.e. more stable and better convergence)
122      !!
123      !! References : Hunke and Dukowicz, JPO97
124      !!              Bouillon et al., Ocean Modelling 2009
125      !!              Bouillon et al., Ocean Modelling 2013
126      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
127      !!-------------------------------------------------------------------
128      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
129      INTEGER                 , INTENT(in   ) ::   Kmm                                   ! ocean time level index
130      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
131      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
132      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   paniso_11 , paniso_12                 ! structure tensor components
133      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   prdg_conv                             ! for ridging
134      !!
135      INTEGER ::   ji, jj       ! dummy loop indices
136      INTEGER ::   jter         ! local integers
137      !
138      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
139      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
140      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
141      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
142      REAl(wp) ::   zbetau, zbetav
143      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
144      REAL(wp) ::   zds2, zdt, zdt2, zdiv, zdiv2, zdsT                  ! temporary scalars
145      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
146      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
147      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
148      !
149      REAL(wp) ::   zintb, zintn                                        ! dummy argument
150      REAL(wp) ::   zfac_x, zfac_y
151      REAL(wp) ::   zshear, zdum1, zdum2
152      REAL(wp) ::   zstressptmp, zstressmtmp, zstress12tmpF             ! anisotropic stress tensor components
153      REAL(wp) ::   zalphar, zalphas                                    ! for mechanical redistribution
154      REAL(wp) ::   zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A  ! for structure tensor evolution
155      !
156      REAL(wp), DIMENSION(jpi,jpj) ::   zstress12tmp                    ! anisotropic stress tensor component for regridding
157      REAL(wp), DIMENSION(jpi,jpj) ::   zyield11, zyield22, zyield12    ! yield surface tensor for history
158      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points
159      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension
160      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
161      !
162      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
163      REAL(wp), DIMENSION(jpi,jpj) ::   zaU  , zaV                      ! ice fraction on U/V points
164      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
165      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
166      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points
167      !
168      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
169      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
170      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
171      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
172      !                                                                 !    ice bottom surface if ice is embedded
173      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses
174      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
175      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
176      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
177      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points
178      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
179      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
180      !
181      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00, zmsk15
182      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
183      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
184
185      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
186      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small
187      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
188      !! --- check convergence
189      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
190      !! --- diags
191      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
192      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p
193      !! --- SIMIP diags
194      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
195      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
196      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
197      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
198      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
199      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)
200      !!-------------------------------------------------------------------
201
202      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology'
203      !
204      ! for diagnostics and convergence tests
205      DO_2D( 1, 1, 1, 1 )
206         zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
207      END_2D
208      IF( nn_rhg_chkcvg > 0 ) THEN
209         DO_2D( 1, 1, 1, 1 )
210            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
211         END_2D
212      ENDIF
213      !
214      !------------------------------------------------------------------------------!
215      ! 0) mask at F points for the ice
216      !------------------------------------------------------------------------------!
217      IF( kt == nit000 ) THEN
218         ! ocean/land mask
219         ALLOCATE( fimask(jpi,jpj) )
220         IF( rn_ishlat == 0._wp ) THEN
221            DO_2D( 0, 0, 0, 0 )
222               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
223            END_2D
224         ELSE
225            DO_2D( 0, 0, 0, 0 )
226               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
227               ! Lateral boundary conditions on velocity (modify fimask)
228               IF( fimask(ji,jj) == 0._wp ) THEN
229                  fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
230                     &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
231               ENDIF
232            END_2D
233         ENDIF
234         CALL lbc_lnk( 'icedyn_rhg_eap', fimask, 'F', 1.0_wp )
235      ENDIF
236
237      !------------------------------------------------------------------------------!
238      ! 1) define some variables and initialize arrays
239      !------------------------------------------------------------------------------!
240      zrhoco = rho0 * rn_cio
241
242      ! ecc2: square of yield ellipse eccenticrity
243      ecc2    = rn_ecc * rn_ecc
244      z1_ecc2 = 1._wp / ecc2
245
246      ! alpha parameters (Bouillon 2009)
247      IF( .NOT. ln_aEVP ) THEN
248         zdtevp   = rDt_ice / REAL( nn_nevp )
249         zalph1 =   2._wp * rn_relast * REAL( nn_nevp )
250         zalph2 = zalph1 * z1_ecc2
251
252         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
253         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
254      ELSE
255         zdtevp   = rdt_ice
256         ! zalpha parameters set later on adaptatively
257      ENDIF
258      z1_dtevp = 1._wp / zdtevp
259
260      ! Initialise stress tensor
261      zs1 (:,:) = pstress1_i (:,:)
262      zs2 (:,:) = pstress2_i (:,:)
263      zs12(:,:) = pstress12_i(:,:)
264
265      ! constants for structure tensor
266      z1_dtevp_A = z1_dtevp/10.0_wp
267      z1dtevpkth = 1._wp / (z1_dtevp_A + 0.00002_wp)
268      zp5kth = 0.5_wp * 0.00002_wp
269
270      ! Ice strength
271      CALL ice_strength
272
273      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
274      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile
275      ELSE                         ;   zkt = 0._wp
276      ENDIF
277      !
278      !------------------------------------------------------------------------------!
279      ! 2) Wind / ocean stress, mass terms, coriolis terms
280      !------------------------------------------------------------------------------!
281      ! sea surface height
282      !    embedded sea ice: compute representative ice top surface
283      !    non-embedded sea ice: use ocean surface for slope calculation
284      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
285
286      DO_2D( 0, 0, 0, 0 )
287
288         ! ice fraction at U-V points
289         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)
290         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)
291
292         ! Ice/snow mass at U-V points
293         zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
294         zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
295         zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
296         zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
297         zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
298
299         ! Ocean currents at U-V points
300         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)
301         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)
302
303         ! Coriolis at T points (m*f)
304         zmf(ji,jj)      = zm1 * ff_t(ji,jj)
305
306         ! dt/m at T points (for alpha and beta coefficients)
307         zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
308
309         ! m/dt
310         zmU_t(ji,jj)    = zmassU * z1_dtevp
311         zmV_t(ji,jj)    = zmassV * z1_dtevp
312
313         ! Drag ice-atm.
314         ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
315         ztauy_ai(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
316
317         ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
318         zspgU(ji,jj)    = - zmassU * grav * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
319         zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
320
321         ! masks
322         zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
323         zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
324
325         ! switches
326         IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
327         ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
328         IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
329         ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
330
331      END_2D
332      CALL lbc_lnk( 'icedyn_rhg_eap', zmf, 'T', 1.0_wp, zdt_m, 'T', 1.0_wp )
333      !
334      !                                  !== Landfast ice parameterization ==!
335      !
336      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
337         DO_2D( 0, 0, 0, 0 )
338            ! ice thickness at U-V points
339            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)
340            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)
341            ! ice-bottom stress at U points
342            zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
343            ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
344            ! ice-bottom stress at V points
345            zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
346            ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
347            ! ice_bottom stress at T points
348            zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) )    ! if grounded icebergs are read: ocean depth = 0
349            tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
350         END_2D
351         CALL lbc_lnk( 'icedyn_rhg_eap', tau_icebfr(:,:), 'T', 1.0_wp )
352         !
353      ELSE                               !-- no landfast
354         DO_2D( 0, 0, 0, 0 )
355            ztaux_base(ji,jj) = 0._wp
356            ztauy_base(ji,jj) = 0._wp
357         END_2D
358      ENDIF
359
360      !------------------------------------------------------------------------------!
361      ! 3) Solution of the momentum equation, iterative procedure
362      !------------------------------------------------------------------------------!
363      !
364      !                                               ! ==================== !
365      DO jter = 1 , nn_nevp                           !    loop over jter    !
366         !                                            ! ==================== !
367         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
368         !
369         ! convergence test
370         IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2  ) THEN
371            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
372               zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
373               zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
374            END_2D
375         ENDIF
376
377         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
378         DO_2D( 1, 0, 1, 0 )
379
380            ! shear at F points
381            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)   &
382               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
383               &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
384
385         END_2D
386
387         DO_2D( 0, 0, 0, 0 )
388
389            ! shear**2 at T points (doc eq. A16)
390            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
391               &   + 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)  &
392               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
393
394            ! divergence at T points
395            zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
396               &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
397               &    ) * r1_e1e2t(ji,jj)
398            zdiv2 = zdiv * zdiv
399
400            ! tension at T points
401            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)   &
402               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
403               &   ) * r1_e1e2t(ji,jj)
404            zdt2 = zdt * zdt
405
406            ! delta at T points
407            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
408
409         END_2D
410         CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1.0_wp )
411
412         ! P/delta at T points
413         DO_2D( 1, 1, 1, 1 )
414            zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
415         END_2D
416
417         DO_2D( 0, 1, 0, 1 )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
418
419             ! shear at T points
420            zdsT = ( zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
421               &   + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
422               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
423
424           ! divergence at T points (duplication to avoid communications)
425            zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
426               &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
427               &    ) * r1_e1e2t(ji,jj)
428
429            ! tension at T points (duplication to avoid communications)
430            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)   &
431               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
432               &   ) * r1_e1e2t(ji,jj)
433
434            ! --- anisotropic stress calculation --- !
435            CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, zdsT, paniso_11(ji,jj), paniso_12(ji,jj), &
436                                    zstressptmp, zstressmtmp, zstress12tmp(ji,jj), strength(ji,jj), zalphar, zalphas)
437
438            ! structure tensor update
439               CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), paniso_11(ji,jj), paniso_12(ji,jj), zmresult11,  zmresult12)
440
441               paniso_11(ji,jj) = (paniso_11(ji,jj)  + 0.5*2.e-5*zdtevp + zmresult11*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit
442               paniso_12(ji,jj) = (paniso_12(ji,jj)                     + zmresult12*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit
443
444            IF (jter == nn_nevp) THEN
445               zyield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp)
446               zyield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp)
447               zyield12(ji,jj) = zstress12tmp(ji,jj)
448               prdg_conv(ji,jj) = -min( zalphar, 0._wp)
449            ENDIF
450
451            ! alpha for aEVP
452            !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
453            !   alpha = beta = sqrt(4*gamma)
454            IF( ln_aEVP ) THEN
455               zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
456               z1_alph1 = 1._wp / ( zalph1 + 1._wp )
457               zalph2   = zalph1
458               z1_alph2 = z1_alph1
459               ! explicit:
460               ! z1_alph1 = 1._wp / zalph1
461               ! z1_alph2 = 1._wp / zalph1
462               ! zalph1 = zalph1 - 1._wp
463               ! zalph2 = zalph1
464            ENDIF
465
466            ! stress at T points (zkt/=0 if landfast)
467            zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zstressptmp ) * z1_alph1
468            zs2(ji,jj) = ( zs2(ji,jj) * zalph1 + zstressmtmp ) * z1_alph1
469         END_2D
470         CALL lbc_lnk( 'icedyn_rhg_eap', zstress12tmp, 'T', 1.0_wp , paniso_11, 'T', 1.0_wp , paniso_12, 'T', 1.0_wp)
471
472        ! Save beta at T-points for further computations
473         IF( ln_aEVP ) THEN
474            DO_2D( 1, 1, 1, 1 )
475               zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
476            END_2D
477         ENDIF
478
479         DO_2D( 1, 0, 1, 0 )
480            ! stress12tmp at F points
481            zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  &
482               &            + zstress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + zstress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  &
483               &            ) * 0.25_wp * r1_e1e2f(ji,jj)
484
485            ! alpha for aEVP
486            IF( ln_aEVP ) THEN
487               zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) )
488               z1_alph2 = 1._wp / ( zalph2 + 1._wp )
489               ! explicit:
490               ! z1_alph2 = 1._wp / zalph2
491               ! zalph2 = zalph2 - 1._wp
492            ENDIF
493
494            ! stress at F points (zkt/=0 if landfast)
495            zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1
496
497         END_2D
498         CALL lbc_lnk( 'icedyn_rhg_eap', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )
499
500         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
501         DO_2D( 0, 0, 0, 0 )
502            !                   !--- U points
503            zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
504               &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
505               &                    ) * r1_e2u(ji,jj)                                                                      &
506               &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
507               &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
508               &                  ) * r1_e1e2u(ji,jj)
509            !
510            !                !--- V points
511            zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
512               &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
513               &                    ) * r1_e1v(ji,jj)                                                                      &
514               &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
515               &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
516               &                  ) * r1_e1e2v(ji,jj)
517            !
518            !                !--- ice currents at U-V point
519            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)
520            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)
521            !
522         END_2D
523         !
524         ! --- Computation of ice velocity --- !
525         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
526         !  Bouillon et al. 2009 (eq 34-35) => stable
527         IF( MOD(jter,2) == 0 ) THEN ! even iterations
528            !
529            DO_2D( 0, 0, 0, 0 )
530               !                 !--- tau_io/(v_oce - v_ice)
531               zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
532                  &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
533               !                 !--- Ocean-to-Ice stress
534               ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
535               !
536               !                 !--- tau_bottom/v_ice
537               zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
538               zTauB = ztauy_base(ji,jj) / zvel
539               !                 !--- OceanBottom-to-Ice stress
540               ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
541               !
542               !                 !--- Coriolis at V-points (energy conserving formulation)
543               zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
544                  &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
545                  &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
546               !
547               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
548               zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
549               !
550               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
551               !                                         1 = sliding friction : TauB < RHS
552               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
553               !
554               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
555                  zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
556                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
557                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
558                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
559                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
560                     &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
561                     &                                    ) / ( zbetav + 1._wp )                                              &
562                     &             ) * 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
563                     &           )   * zmsk00y(ji,jj)
564               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
565                  v_ice(ji,jj) = ( (         rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
566                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
567                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
568                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
569                     &              ) * 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
570                     &            )   * zmsk00y(ji,jj)
571               ENDIF
572            END_2D
573            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
574            !
575#if defined key_agrif
576!!          CALL agrif_interp_ice( 'V', jter, nn_nevp )
577            CALL agrif_interp_ice( 'V' )
578#endif
579            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
580            !
581            DO_2D( 0, 0, 0, 0 )
582               !                 !--- tau_io/(u_oce - u_ice)
583               zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
584                  &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
585               !                 !--- Ocean-to-Ice stress
586               ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
587               !
588               !                 !--- tau_bottom/u_ice
589               zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
590               zTauB = ztaux_base(ji,jj) / zvel
591               !                 !--- OceanBottom-to-Ice stress
592               ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
593               !
594               !                 !--- Coriolis at U-points (energy conserving formulation)
595               zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
596                  &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
597                  &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
598               !
599               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
600               zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
601               !
602               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
603               !                                         1 = sliding friction : TauB < RHS
604               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
605               !
606               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
607                  zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
608                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
609                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
610                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
611                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
612                     &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
613                     &                                    ) / ( zbetau + 1._wp )                                              &
614                     &             ) * 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
615                     &           )   * zmsk00x(ji,jj)
616               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
617                  u_ice(ji,jj) = ( (         rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
618                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
619                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
620                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
621                     &              ) * 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
622                     &            )   * zmsk00x(ji,jj)
623               ENDIF
624            END_2D
625            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
626            !
627#if defined key_agrif
628!!          CALL agrif_interp_ice( 'U', jter, nn_nevp )
629            CALL agrif_interp_ice( 'U' )
630#endif
631            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
632            !
633         ELSE ! odd iterations
634            !
635            DO_2D( 0, 0, 0, 0 )
636               !                 !--- tau_io/(u_oce - u_ice)
637               zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
638                  &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
639               !                 !--- Ocean-to-Ice stress
640               ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
641               !
642               !                 !--- tau_bottom/u_ice
643               zvel  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
644               zTauB = ztaux_base(ji,jj) / zvel
645               !                 !--- OceanBottom-to-Ice stress
646               ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
647               !
648               !                 !--- Coriolis at U-points (energy conserving formulation)
649               zCorU(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
650                  &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
651                  &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
652               !
653               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
654               zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
655               !
656               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
657               !                                         1 = sliding friction : TauB < RHS
658               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
659               !
660               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
661                  zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
662                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
663                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
664                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
665                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
666                     &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
667                     &                                    ) / ( zbetau + 1._wp )                                              &
668                     &             ) * 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
669                     &           )   * zmsk00x(ji,jj)
670               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
671                  u_ice(ji,jj) = ( (         rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
672                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
673                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
674                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
675                     &              ) * 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
676                     &            )   * zmsk00x(ji,jj)
677               ENDIF
678            END_2D
679            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1.0_wp )
680            !
681#if defined key_agrif
682!!          CALL agrif_interp_ice( 'U', jter, nn_nevp )
683            CALL agrif_interp_ice( 'U' )
684#endif
685            IF( ln_bdy )   CALL bdy_ice_dyn( 'U' )
686            !
687            DO_2D( 0, 0, 0, 0 )
688               !                 !--- tau_io/(v_oce - v_ice)
689               zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
690                  &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
691               !                 !--- Ocean-to-Ice stress
692               ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
693               !
694               !                 !--- tau_bottom/v_ice
695               zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
696               zTauB = ztauy_base(ji,jj) / zvel
697               !                 !--- OceanBottom-to-Ice stress
698               ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
699               !
700               !                 !--- Coriolis at v-points (energy conserving formulation)
701               zCorV(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
702                  &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
703                  &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
704               !
705               !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
706               zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
707               !
708               !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
709               !                                         1 = sliding friction : TauB < RHS
710               rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
711               !
712               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
713                  zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
714                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
715                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
716                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
717                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
718                     &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
719                     &                                    ) / ( zbetav + 1._wp )                                              &
720                     &             ) * 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
721                     &           )   * zmsk00y(ji,jj)
722               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
723                  v_ice(ji,jj) = ( (         rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
724                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
725                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast
726                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0
727                     &              ) * 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
728                     &            )   * zmsk00y(ji,jj)
729               ENDIF
730            END_2D
731            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1.0_wp )
732            !
733#if defined key_agrif
734!!          CALL agrif_interp_ice( 'V', jter, nn_nevp )
735            CALL agrif_interp_ice( 'V' )
736#endif
737            IF( ln_bdy )   CALL bdy_ice_dyn( 'V' )
738            !
739         ENDIF
740
741         ! convergence test
742         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 )
743         !
744         !                                                ! ==================== !
745      END DO                                              !  end loop over jter  !
746      !                                                   ! ==================== !
747      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta )
748      !
749      CALL lbc_lnk( 'icedyn_rhg_eap', prdg_conv, 'T', 1.0_wp )  ! only need this in ridging module after rheology completed
750      !
751      !------------------------------------------------------------------------------!
752      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
753      !------------------------------------------------------------------------------!
754      DO_2D( 1, 0, 1, 0 )
755
756         ! shear at F points
757         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)   &
758            &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
759            &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
760
761      END_2D
762
763      DO_2D( 0, 0, 0, 0 )
764
765         ! tension**2 at T points
766         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)   &
767            &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
768            &   ) * r1_e1e2t(ji,jj)
769         zdt2 = zdt * zdt
770
771         zten_i(ji,jj) = zdt
772
773         ! shear**2 at T points (doc eq. A16)
774         zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
775            &   + 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)  &
776            &   ) * 0.25_wp * r1_e1e2t(ji,jj)
777
778         ! shear at T points
779         pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
780
781         ! divergence at T points
782         pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
783            &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
784            &             ) * r1_e1e2t(ji,jj)
785
786         ! delta at T points
787         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
788         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
789         pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
790
791      END_2D
792      CALL lbc_lnk( 'icedyn_rhg_eap', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp, &
793         &                              zten_i, 'T', 1.0_wp, zs1    , 'T', 1.0_wp, zs2     , 'T', 1.0_wp, &
794         &                                zs12, 'F', 1.0_wp )
795
796      ! --- Store the stress tensor for the next time step --- !
797      pstress1_i (:,:) = zs1 (:,:)
798      pstress2_i (:,:) = zs2 (:,:)
799      pstress12_i(:,:) = zs12(:,:)
800      !
801
802      !------------------------------------------------------------------------------!
803      ! 5) diagnostics
804      !------------------------------------------------------------------------------!
805      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
806      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
807         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
808         !
809         CALL lbc_lnk( 'icedyn_rhg_eap', ztaux_oi, 'U', -1.0_wp, ztauy_oi, 'V', -1.0_wp, ztaux_ai, 'U', -1.0_wp, &
810            &                            ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
811         !
812         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
813         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
814         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
815         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
816         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
817         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
818      ENDIF
819
820      ! --- divergence, shear and strength --- !
821      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
822      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
823      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
824      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
825
826      ! --- Stress tensor invariants (SIMIP diags) --- !
827      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
828         !
829         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
830         !
831         DO_2D( 1, 1, 1, 1 )
832
833            ! Ice stresses
834            ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
835            ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
836            ! I know, this can be confusing...
837            zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl )
838            zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
839            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
840            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
841
842            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
843            zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure
844            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
845
846         END_2D
847         !
848         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
849         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
850         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
851
852         DEALLOCATE ( zsig_I, zsig_II )
853
854      ENDIF
855
856      ! --- Normalized stress tensor principal components --- !
857      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
858      ! Recommendation 1 : we use ice strength, not replacement pressure
859      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
860      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
861         !
862         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
863         !
864         DO_2D( 1, 1, 1, 1 )
865
866            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
867            !                        and **deformations** at current iterates
868            !                        following Lemieux & Dupont (2020)
869            zfac             =   zp_delt(ji,jj)
870            zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
871            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
872            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
873
874            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
875            zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure
876            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
877
878            ! Normalized  principal stresses (used to display the ellipse)
879            z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )
880            zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
881            zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
882         END_2D
883         !
884         CALL iom_put( 'sig1_pnorm' , zsig1_p )
885         CALL iom_put( 'sig2_pnorm' , zsig2_p )
886
887         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
888
889      ENDIF
890
891      ! --- yieldcurve --- !
892      IF( iom_use('yield11') .OR. iom_use('yield12') .OR. iom_use('yield22')) THEN
893
894         CALL lbc_lnk( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp )
895
896         CALL iom_put( 'yield11', zyield11 * zmsk00 )
897         CALL iom_put( 'yield22', zyield22 * zmsk00 )
898         CALL iom_put( 'yield12', zyield12 * zmsk00 )
899      ENDIF
900
901      ! --- anisotropy tensor --- !
902      IF( iom_use('aniso') ) THEN
903         CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp )
904         CALL iom_put( 'aniso' , paniso_11 * zmsk00 )
905      ENDIF
906
907      ! --- SIMIP --- !
908      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
909         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
910         !
911         CALL lbc_lnk( 'icedyn_rhg_eap', zspgU, 'U', -1.0_wp, zspgV, 'V', -1.0_wp, &
912            &                            zCorU, 'U', -1.0_wp, zCorV, 'V', -1.0_wp, &
913            &                              zfU, 'U', -1.0_wp,   zfV, 'V', -1.0_wp )
914
915         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
916         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
917         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
918         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
919         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
920         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
921      ENDIF
922
923      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
924         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
925         !
926         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
927            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
928         !
929         DO_2D( 0, 0, 0, 0 )
930            ! 2D ice mass, snow mass, area transport arrays (X, Y)
931            zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
932            zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
933
934            zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
935            zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
936
937            zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
938            zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
939
940            zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
941            zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
942
943         END_2D
944
945         CALL lbc_lnk( 'icedyn_rhg_eap', zdiag_xmtrp_ice, 'U', -1.0_wp, zdiag_ymtrp_ice, 'V', -1.0_wp, &
946            &                            zdiag_xmtrp_snw, 'U', -1.0_wp, zdiag_ymtrp_snw, 'V', -1.0_wp, &
947            &                            zdiag_xatrp    , 'U', -1.0_wp, zdiag_yatrp    , 'V', -1.0_wp )
948
949         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
950         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport
951         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
952         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
953         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
954         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
955
956         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
957            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
958
959      ENDIF
960      !
961      ! --- convergence tests --- !
962      IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
963         IF( iom_use('uice_cvg') ) THEN
964            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
965               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
966                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
967            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
968               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
969                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
970            ENDIF
971         ENDIF
972      ENDIF
973      !
974   END SUBROUTINE ice_dyn_rhg_eap
975
976
977   SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 )
978      !!----------------------------------------------------------------------
979      !!                    ***  ROUTINE rhg_cvg_eap  ***
980      !!
981      !! ** Purpose :   check convergence of oce rheology
982      !!
983      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity
984      !!                during the sub timestepping of rheology so as:
985      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) )
986      !!                This routine is called every sub-iteration, so it is cpu expensive
987      !!
988      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)
989      !!----------------------------------------------------------------------
990      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
991      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
992      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pmsk15
993      !!
994      INTEGER           ::   it, idtime, istatus
995      INTEGER           ::   ji, jj          ! dummy loop indices
996      REAL(wp)          ::   zresm           ! local real
997      CHARACTER(len=20) ::   clname
998      !!----------------------------------------------------------------------
999
1000      ! create file
1001      IF( kt == nit000 .AND. kiter == 1 ) THEN
1002         !
1003         IF( lwp ) THEN
1004            WRITE(numout,*)
1005            WRITE(numout,*) 'rhg_cvg_eap : ice rheology convergence control'
1006            WRITE(numout,*) '~~~~~~~'
1007         ENDIF
1008         !
1009         IF( lwm ) THEN
1010            clname = 'ice_cvg.nc'
1011            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
1012            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
1013            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
1014            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid )
1015            istatus = NF90_ENDDEF(ncvgid)
1016         ENDIF
1017         !
1018      ENDIF
1019
1020      ! time
1021      it = ( kt - 1 ) * kitermax + kiter
1022
1023      ! convergence
1024      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
1025         zresm = 0._wp
1026      ELSE
1027         zresm = 0._wp
1028         DO_2D( 0, 0, 0, 0 )
1029            zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
1030               &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
1031         END_2D
1032         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
1033      ENDIF
1034
1035      IF( lwm ) THEN
1036         ! write variables
1037         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
1038         ! close file
1039         IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid)
1040      ENDIF
1041
1042   END SUBROUTINE rhg_cvg_eap
1043
1044
1045   SUBROUTINE update_stress_rdg( ksub, kndte, pdivu, ptension, pshear, pa11, pa12, &
1046      &                          pstressp,  pstressm, pstress12, pstrength, palphar, palphas )
1047      !!---------------------------------------------------------------------
1048      !!                   ***  SUBROUTINE update_stress_rdg  ***
1049      !!
1050      !! ** Purpose :   Updates the stress depending on values of strain rate and structure
1051      !!                tensor and for the last subcycle step it computes closing and sliding rate
1052      !!---------------------------------------------------------------------
1053      INTEGER,  INTENT(in   ) ::   ksub, kndte
1054      REAL(wp), INTENT(in   ) ::   pstrength
1055      REAL(wp), INTENT(in   ) ::   pdivu, ptension, pshear
1056      REAL(wp), INTENT(in   ) ::   pa11, pa12
1057      REAL(wp), INTENT(  out) ::   pstressp, pstressm, pstress12
1058      REAL(wp), INTENT(  out) ::   palphar, palphas
1059
1060      INTEGER  ::   kx ,ky, ka
1061
1062      REAL(wp) ::   zstemp11r, zstemp12r, zstemp22r
1063      REAL(wp) ::   zstemp11s, zstemp12s, zstemp22s
1064      REAL(wp) ::   za22, zQd11Qd11, zQd11Qd12, zQd12Qd12
1065      REAL(wp) ::   zQ11Q11, zQ11Q12, zQ12Q12
1066      REAL(wp) ::   zdtemp11, zdtemp12, zdtemp22
1067      REAL(wp) ::   zrotstemp11r, zrotstemp12r, zrotstemp22r
1068      REAL(wp) ::   zrotstemp11s, zrotstemp12s, zrotstemp22s
1069      REAL(wp) ::   zsig11, zsig12, zsig22
1070      REAL(wp) ::   zsgprm11, zsgprm12, zsgprm22
1071      REAL(wp) ::   zAngle_denom_gamma,  zAngle_denom_alpha
1072      REAL(wp) ::   zTany_1, zTany_2
1073      REAL(wp) ::   zx, zy, zkxw, kyw, kaw
1074      REAL(wp) ::   zinvdx, zinvdy, zinvda
1075      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime
1076
1077      REAL(wp), PARAMETER ::   ppkfriction = 0.45_wp
1078      ! Factor to maintain the same stress as in EVP (see Section 3)
1079      ! Can be set to 1 otherwise
1080!     REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp/(1._wp+ppkfriction*ppkfriction)
1081      REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp
1082
1083      ! next statement uses pphi set in main module (icedyn_rhg_eap)
1084      REAL(wp), PARAMETER ::   ppinvsin = 1._wp/sin(2._wp*pphi) * ppinvstressconviso
1085
1086      ! compute eigenvalues, eigenvectors and angles for structure tensor, strain
1087      ! rates
1088
1089      ! 1) structure tensor
1090      za22 = 1._wp-pa11
1091      zQ11Q11 = 1._wp
1092      zQ12Q12 = rsmall
1093      zQ11Q12 = rsmall
1094
1095      ! gamma: angle between general coordiantes and principal axis of A
1096      ! here Tan2gamma = 2 a12 / (a11 - a22)
1097      IF((ABS(pa11 - za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN
1098         zAngle_denom_gamma = 1._wp/sqrt( ( pa11 - za22 )*( pa11 - za22) + &
1099                              4._wp*pa12*pa12 )
1100
1101         zQ11Q11 = 0.5_wp + ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Cos^2
1102         zQ12Q12 = 0.5_wp - ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Sin^2
1103         zQ11Q12 = pa12*zAngle_denom_gamma                     !CosSin
1104      ENDIF
1105
1106      ! rotation Q*atemp*Q^T
1107      zatempprime = zQ11Q11*pa11 + 2.0_wp*zQ11Q12*pa12 + zQ12Q12*za22
1108
1109      ! make first principal value the largest
1110      zatempprime = max(zatempprime, 1.0_wp - zatempprime)
1111
1112      ! 2) strain rate
1113      zdtemp11 = 0.5_wp*(pdivu + ptension)
1114      zdtemp12 = pshear*0.5_wp
1115      zdtemp22 = 0.5_wp*(pdivu - ptension)
1116
1117      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22)
1118
1119      zQd11Qd11 = 1.0_wp
1120      zQd12Qd12 = rsmall
1121      zQd11Qd12 = rsmall
1122
1123      IF((ABS( zdtemp11 - zdtemp22) > rsmall).OR. (ABS(zdtemp12) > rsmall)) THEN
1124
1125         zAngle_denom_alpha = 1.0_wp/sqrt( ( zdtemp11 - zdtemp22 )* &
1126                             ( zdtemp11 - zdtemp22 ) + 4.0_wp*zdtemp12*zdtemp12)
1127
1128         zQd11Qd11 = 0.5_wp + ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Cos^2
1129         zQd12Qd12 = 0.5_wp - ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Sin^2
1130         zQd11Qd12 = zdtemp12*zAngle_denom_alpha !CosSin
1131      ENDIF
1132
1133      zdtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22
1134      zdtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22
1135      ! In cos and sin values
1136      zx = 0._wp
1137      IF ((ABS(zdtemp1) > rsmall).OR.(ABS(zdtemp2) > rsmall)) THEN
1138         zx = atan2(zdtemp2,zdtemp1)
1139      ENDIF
1140
1141      ! to ensure the angle lies between pi/4 and 9 pi/4
1142      IF (zx < rpi*0.25_wp) zx = zx + rpi*2.0_wp
1143
1144      ! y: angle between major principal axis of strain rate and structure
1145      ! tensor
1146      ! y = gamma - alpha
1147      ! Expressesed componently with
1148      ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha)
1149
1150      zTany_1 = zQ11Q12 - zQd11Qd12
1151      zTany_2 = zQ11Q11 - zQd12Qd12
1152
1153      zy = 0._wp
1154
1155      IF ((ABS(zTany_1) > rsmall).OR.(ABS(zTany_2) > rsmall)) THEN
1156         zy = atan2(zTany_1,zTany_2)
1157      ENDIF
1158
1159      ! to make sure y is between 0 and pi
1160      IF (zy > rpi) zy = zy - rpi
1161      IF (zy < 0)  zy = zy + rpi
1162
1163      ! 3) update anisotropic stress tensor
1164      zinvdx = real(nx_yield-1,kind=wp)/rpi
1165      zinvdy = real(ny_yield-1,kind=wp)/rpi
1166      zinvda = 2._wp*real(na_yield-1,kind=wp)
1167
1168      ! % need 8 coords and 8 weights
1169      ! % range in kx
1170      kx  = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1
1171      !!clem kx  = MAX( 1, MIN( nx_yield-1, INT((zx-rpi*0.25_wp-rpi)*zinvdx) + 1  ) )
1172      zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx
1173
1174      ky  = int(zy*zinvdy) + 1
1175      !!clem ky  = MAX( 1, MIN( ny_yield-1, INT(zy*zinvdy) + 1 ) )
1176      kyw = ky - zy*zinvdy
1177
1178      ka  = int((zatempprime-0.5_wp)*zinvda) + 1
1179      !!clem ka  = MAX( 1, MIN( na_yield-1, INT((zatempprime-0.5_wp)*zinvda) + 1 ) )
1180      kaw = ka - (zatempprime-0.5_wp)*zinvda
1181
1182      ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013)
1183!!$      zstemp11r =  zkxw  * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) &
1184!!$        & + (1._wp-zkxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) &
1185!!$        & + zkxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) &
1186!!$        & + zkxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) &
1187!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) &
1188!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) &
1189!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) &
1190!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1)
1191!!$      zstemp12r =  zkxw  * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) &
1192!!$        & + (1._wp-zkxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) &
1193!!$        & + zkxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) &
1194!!$        & + zkxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) &
1195!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) &
1196!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) &
1197!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) &
1198!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1)
1199!!$      zstemp22r =  zkxw  * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) &
1200!!$        & + (1._wp-zkxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) &
1201!!$        & + zkxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) &
1202!!$        & + zkxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) &
1203!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) &
1204!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) &
1205!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) &
1206!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1)
1207!!$
1208!!$      zstemp11s =  zkxw  * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) &
1209!!$        & + (1._wp-zkxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) &
1210!!$        & + zkxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) &
1211!!$        & + zkxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) &
1212!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) &
1213!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) &
1214!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) &
1215!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1)
1216!!$      zstemp12s =  zkxw  * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) &
1217!!$        & + (1._wp-zkxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) &
1218!!$        & + zkxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) &
1219!!$        & + zkxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) &
1220!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) &
1221!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) &
1222!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) &
1223!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1)
1224!!$      zstemp22s =  zkxw  * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) &
1225!!$        & + (1._wp-zkxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) &
1226!!$        & + zkxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) &
1227!!$        & + zkxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) &
1228!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) &
1229!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) &
1230!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) &
1231!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1)
1232
1233      zstemp11r = s11r(kx,ky,ka)
1234      zstemp12r = s12r(kx,ky,ka)
1235      zstemp22r = s22r(kx,ky,ka)
1236      zstemp11s = s11s(kx,ky,ka)
1237      zstemp12s = s12s(kx,ky,ka)
1238      zstemp22s = s22s(kx,ky,ka)
1239
1240
1241      ! Calculate mean ice stress over a collection of floes (Equation 3 in
1242      ! Tsamados 2013)
1243
1244      zsig11  = pstrength*(zstemp11r + ppkfriction*zstemp11s) * ppinvsin
1245      zsig12  = pstrength*(zstemp12r + ppkfriction*zstemp12s) * ppinvsin
1246      zsig22  = pstrength*(zstemp22r + ppkfriction*zstemp22s) * ppinvsin
1247
1248      ! Back - rotation of the stress from principal axes into general coordinates
1249
1250      ! Update stress
1251      zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 -      2._wp*zQ11Q12 *zsig12
1252      zsgprm12 = zQ11Q12*zsig11 - zQ11Q12*zsig22 + (zQ11Q11 - zQ12Q12)*zsig12
1253      zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 +      2._wp*zQ11Q12 *zsig12
1254
1255      pstressp  = zsgprm11 + zsgprm22
1256      pstress12 = zsgprm12
1257      pstressm  = zsgprm11 - zsgprm22
1258
1259      ! Compute ridging and sliding functions in general coordinates
1260      ! (Equation 11 in Tsamados 2013)
1261      IF (ksub == kndte) THEN
1262         zrotstemp11r = zQ11Q11*zstemp11r - 2._wp*zQ11Q12* zstemp12r &
1263                      + zQ12Q12*zstemp22r
1264         zrotstemp12r = zQ11Q11*zstemp12r +       zQ11Q12*(zstemp11r-zstemp22r) &
1265                      - zQ12Q12*zstemp12r
1266         zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r &
1267                      + zQ11Q11*zstemp22r
1268
1269         zrotstemp11s = zQ11Q11*zstemp11s - 2._wp*zQ11Q12* zstemp12s &
1270                      + zQ12Q12*zstemp22s
1271         zrotstemp12s = zQ11Q11*zstemp12s +       zQ11Q12*(zstemp11s-zstemp22s) &
1272                      - zQ12Q12*zstemp12s
1273         zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s &
1274                      + zQ11Q11*zstemp22s
1275
1276         palphar =  zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 &
1277                  + zrotstemp22r*zdtemp22
1278         palphas =  zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 &
1279                  + zrotstemp22s*zdtemp22
1280      ENDIF
1281   END SUBROUTINE update_stress_rdg
1282
1283!=======================================================================
1284
1285
1286   SUBROUTINE calc_ffrac( pstressp, pstressm, pstress12, pa11, pa12, &
1287      &                   pmresult11, pmresult12 )
1288      !!---------------------------------------------------------------------
1289      !!                   ***  ROUTINE calc_ffrac  ***
1290      !!
1291      !! ** Purpose :   Computes term in evolution equation for structure tensor
1292      !!                which determines the ice floe re-orientation due to fracture
1293      !! ** Method :    Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2
1294      !!---------------------------------------------------------------------
1295      REAL (wp), INTENT(in)  ::   pstressp, pstressm, pstress12, pa11, pa12
1296      REAL (wp), INTENT(out) ::   pmresult11, pmresult12
1297
1298      ! local variables
1299
1300      REAL (wp) ::   zsigma11, zsigma12, zsigma22  ! stress tensor elements
1301      REAL (wp) ::   zAngle_denom        ! angle with principal component axis
1302      REAL (wp) ::   zsigma_1, zsigma_2   ! principal components of stress
1303      REAL (wp) ::   zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12
1304
1305!!$   REAL (wp), PARAMETER ::   ppkfrac = 0.0001_wp   ! rate of fracture formation
1306      REAL (wp), PARAMETER ::   ppkfrac = 1.e-3_wp    ! rate of fracture formation
1307      REAL (wp), PARAMETER ::   ppthreshold = 0.3_wp  ! critical confinement ratio
1308      !!---------------------------------------------------------------
1309      !
1310      zsigma11 = 0.5_wp*(pstressp+pstressm)
1311      zsigma12 = pstress12
1312      zsigma22 = 0.5_wp*(pstressp-pstressm)
1313
1314      ! Here's the change - no longer calculate gamma,
1315      ! use trig formulation, angle signs are kept correct, don't worry
1316
1317      ! rotate tensor to get into sigma principal axis
1318
1319      ! here Tan2gamma = 2 sig12 / (sig11 - sig12)
1320      ! add rsmall to the denominator to stop 1/0 errors, makes very little
1321      ! error to the calculated angles < rsmall
1322
1323      zQ11Q11 = 0.1_wp
1324      zQ12Q12 = rsmall
1325      zQ11Q12 = rsmall
1326
1327      IF((ABS( zsigma11 - zsigma22) > rsmall).OR.(ABS(zsigma12) > rsmall)) THEN
1328
1329         zAngle_denom = 1.0_wp/sqrt( ( zsigma11 - zsigma22 )*( zsigma11 - &
1330                       zsigma22 ) + 4.0_wp*zsigma12*zsigma12)
1331
1332         zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Cos^2
1333         zQ12Q12 = 0.5_wp - ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Sin^2
1334         zQ11Q12 = zsigma12*zAngle_denom                      !CosSin
1335      ENDIF
1336
1337      zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12 + zQ12Q12*zsigma22 ! S(1,1)
1338      zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12 + zQ11Q11*zsigma22 ! S(2,2)
1339
1340      ! Pure divergence
1341      IF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 >= 0.0_wp))  THEN
1342         pmresult11 = 0.0_wp
1343         pmresult12 = 0.0_wp
1344
1345      ! Unconfined compression: cracking of blocks not along the axial splitting
1346      ! direction
1347      ! which leads to the loss of their shape, so we again model it through diffusion
1348      ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp))  THEN
1349         pmresult11 = - ppkfrac * (pa11 - zQ12Q12)
1350         pmresult12 = - ppkfrac * (pa12 + zQ11Q12)
1351
1352      ! Shear faulting
1353      ELSEIF (zsigma_2 == 0.0_wp) THEN
1354         pmresult11 = 0.0_wp
1355         pmresult12 = 0.0_wp
1356      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= ppthreshold)) THEN
1357         pmresult11 = - ppkfrac * (pa11 - zQ12Q12)
1358         pmresult12 = - ppkfrac * (pa12 + zQ11Q12)
1359
1360      ! Horizontal spalling
1361      ELSE
1362         pmresult11 = 0.0_wp
1363         pmresult12 = 0.0_wp
1364      ENDIF
1365
1366   END SUBROUTINE calc_ffrac
1367
1368
1369   SUBROUTINE rhg_eap_rst( cdrw, kt )
1370      !!---------------------------------------------------------------------
1371      !!                   ***  ROUTINE rhg_eap_rst  ***
1372      !!
1373      !! ** Purpose :   Read or write RHG file in restart file
1374      !!
1375      !! ** Method  :   use of IOM library
1376      !!----------------------------------------------------------------------
1377      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1378      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
1379      !
1380      INTEGER  ::   iter                      ! local integer
1381      INTEGER  ::   id1, id2, id3, id4, id5   ! local integers
1382      INTEGER  ::   ix, iy, ip, iz, n, ia     ! local integers
1383
1384      INTEGER, PARAMETER            ::    nz = 100
1385
1386      REAL(wp) ::   ainit, xinit, yinit, pinit, zinit
1387      REAL(wp) ::   da, dx, dy, dp, dz, a1
1388
1389      !!clem
1390      REAL(wp) ::   zw1, zw2, zfac, ztemp
1391      REAL(wp) ::   zidx, zidy, zidz
1392      REAL(wp) ::   zsaak(6)                  ! temporary array
1393
1394      REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp
1395      !!----------------------------------------------------------------------
1396      !
1397      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
1398         !                                   ! ---------------
1399         IF( ln_rstart ) THEN                   !* Read the restart file
1400            !
1401            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
1402            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
1403            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
1404            id4 = iom_varid( numrir, 'aniso_11'  , ldstop = .FALSE. )
1405            id5 = iom_varid( numrir, 'aniso_12'  , ldstop = .FALSE. )
1406            !
1407            IF( MIN( id1, id2, id3, id4, id5 ) > 0 ) THEN      ! fields exist
1408               CALL iom_get( numrir, jpdom_auto, 'stress1_i' , stress1_i , cd_type = 'T' )
1409               CALL iom_get( numrir, jpdom_auto, 'stress2_i' , stress2_i , cd_type = 'T' )
1410               CALL iom_get( numrir, jpdom_auto, 'stress12_i', stress12_i, cd_type = 'F' )
1411               CALL iom_get( numrir, jpdom_auto, 'aniso_11'  , aniso_11  , cd_type = 'T' )
1412               CALL iom_get( numrir, jpdom_auto, 'aniso_12'  , aniso_12  , cd_type = 'T' )
1413            ELSE                                     ! start rheology from rest
1414               IF(lwp) WRITE(numout,*)
1415               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
1416               stress1_i (:,:) = 0._wp
1417               stress2_i (:,:) = 0._wp
1418               stress12_i(:,:) = 0._wp
1419               aniso_11  (:,:) = 0.5_wp
1420               aniso_12  (:,:) = 0._wp
1421            ENDIF
1422         ELSE                                   !* Start from rest
1423            IF(lwp) WRITE(numout,*)
1424            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
1425            stress1_i (:,:) = 0._wp
1426            stress2_i (:,:) = 0._wp
1427            stress12_i(:,:) = 0._wp
1428            aniso_11  (:,:) = 0.5_wp
1429            aniso_12  (:,:) = 0._wp
1430         ENDIF
1431         !
1432
1433         da = 0.5_wp/real(na_yield-1,kind=wp)
1434         ainit = 0.5_wp - da
1435         dx = rpi/real(nx_yield-1,kind=wp)
1436         xinit = rpi + 0.25_wp*rpi - dx
1437         dz = rpi/real(nz,kind=wp)
1438         zinit = -rpi*0.5_wp
1439         dy = rpi/real(ny_yield-1,kind=wp)
1440         yinit = -dy
1441
1442         s11r(:,:,:) = 0._wp
1443         s12r(:,:,:) = 0._wp
1444         s22r(:,:,:) = 0._wp
1445         s11s(:,:,:) = 0._wp
1446         s12s(:,:,:) = 0._wp
1447         s22s(:,:,:) = 0._wp
1448
1449!!$         DO ia=1,na_yield
1450!!$            DO ix=1,nx_yield
1451!!$               DO iy=1,ny_yield
1452!!$                  s11r(ix,iy,ia) = 0._wp
1453!!$                  s12r(ix,iy,ia) = 0._wp
1454!!$                  s22r(ix,iy,ia) = 0._wp
1455!!$                  s11s(ix,iy,ia) = 0._wp
1456!!$                  s12s(ix,iy,ia) = 0._wp
1457!!$                  s22s(ix,iy,ia) = 0._wp
1458!!$                  IF (ia <= na_yield-1) THEN
1459!!$                     DO iz=1,nz
1460!!$                        s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
1461!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
1462!!$                           s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
1463!!$                        s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
1464!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
1465!!$                           s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
1466!!$                        s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
1467!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
1468!!$                           s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
1469!!$                        s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
1470!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
1471!!$                           s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
1472!!$                        s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
1473!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
1474!!$                           s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
1475!!$                        s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
1476!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
1477!!$                           s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
1478!!$                     ENDDO
1479!!$                     IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp
1480!!$                     IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp
1481!!$                     IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp
1482!!$                     IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp
1483!!$                     IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp
1484!!$                     IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp
1485!!$                  ELSE
1486!!$                     s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
1487!!$                     s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
1488!!$                     s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
1489!!$                     s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
1490!!$                     s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
1491!!$                     s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
1492!!$                     IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp
1493!!$                     IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp
1494!!$                     IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp
1495!!$                     IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp
1496!!$                     IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp
1497!!$                     IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp
1498!!$                  ENDIF
1499!!$               ENDDO
1500!!$            ENDDO
1501!!$         ENDDO
1502
1503         !! faster but still very slow => to be improved
1504         zfac = dz/sin(2._wp*pphi)
1505         DO ia = 1, na_yield-1
1506            zw1 = w1(ainit+ia*da)
1507            zw2 = w2(ainit+ia*da)
1508            DO iz = 1, nz
1509               zidz = zinit+iz*dz
1510               ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz))
1511               DO iy = 1, ny_yield
1512                  zidy = yinit+iy*dy
1513                  DO ix = 1, nx_yield
1514                     zidx = xinit+ix*dx
1515                     call all_skr_sks(zidx,zidy,zidz,zsaak)
1516                     zsaak = ztemp*zsaak*zfac
1517                     s11r(ix,iy,ia) = s11r(ix,iy,ia) + zsaak(1)
1518                     s12r(ix,iy,ia) = s12r(ix,iy,ia) + zsaak(2)
1519                     s22r(ix,iy,ia) = s22r(ix,iy,ia) + zsaak(3)
1520                     s11s(ix,iy,ia) = s11s(ix,iy,ia) + zsaak(4)
1521                     s12s(ix,iy,ia) = s12s(ix,iy,ia) + zsaak(5)
1522                     s22s(ix,iy,ia) = s22s(ix,iy,ia) + zsaak(6)
1523                  END DO
1524               END DO
1525            END DO
1526         END DO
1527         zfac = 1._wp/sin(2._wp*pphi)
1528         ia = na_yield
1529         DO iy = 1, ny_yield
1530            zidy = yinit+iy*dy
1531            DO ix = 1, nx_yield
1532               zidx = xinit+ix*dx
1533               call all_skr_sks(zidx,zidy,0._wp,zsaak)
1534               zsaak = 0.5_wp*zsaak*zfac
1535               s11r(ix,iy,ia) = zsaak(1)
1536               s12r(ix,iy,ia) = zsaak(2)
1537               s22r(ix,iy,ia) = zsaak(3)
1538               s11s(ix,iy,ia) = zsaak(4)
1539               s12s(ix,iy,ia) = zsaak(5)
1540               s22s(ix,iy,ia) = zsaak(6)
1541            ENDDO
1542         ENDDO
1543         WHERE (ABS(s11r(:,:,:)) < eps6) s11r(:,:,:) = 0._wp
1544         WHERE (ABS(s12r(:,:,:)) < eps6) s12r(:,:,:) = 0._wp
1545         WHERE (ABS(s22r(:,:,:)) < eps6) s22r(:,:,:) = 0._wp
1546         WHERE (ABS(s11s(:,:,:)) < eps6) s11s(:,:,:) = 0._wp
1547         WHERE (ABS(s12s(:,:,:)) < eps6) s12s(:,:,:) = 0._wp
1548         WHERE (ABS(s22s(:,:,:)) < eps6) s22s(:,:,:) = 0._wp
1549
1550
1551      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1552         !                                   ! -------------------
1553         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
1554         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
1555         !
1556         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
1557         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
1558         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
1559         CALL iom_rstput( iter, nitrst, numriw, 'aniso_11'  , aniso_11   )
1560         CALL iom_rstput( iter, nitrst, numriw, 'aniso_12'  , aniso_12   )
1561         !
1562      ENDIF
1563      !
1564   END SUBROUTINE rhg_eap_rst
1565
1566
1567   FUNCTION w1(pa)
1568      !!-------------------------------------------------------------------
1569      !! Function : w1 (see Gaussian function psi in Tsamados et al 2013)
1570      !!-------------------------------------------------------------------
1571      REAL(wp), INTENT(in   ) ::   pa
1572      REAL(wp) ::   w1
1573      !!-------------------------------------------------------------------
1574
1575      w1 = -   223.87569446_wp &
1576       &   +  2361.21986630_wp*pa &
1577       &   - 10606.56079975_wp*pa*pa &
1578       &   + 26315.50025642_wp*pa*pa*pa &
1579       &   - 38948.30444297_wp*pa*pa*pa*pa &
1580       &   + 34397.72407466_wp*pa*pa*pa*pa*pa &
1581       &   - 16789.98003081_wp*pa*pa*pa*pa*pa*pa &
1582       &   +  3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa
1583
1584   END FUNCTION w1
1585
1586
1587   FUNCTION w2(pa)
1588      !!-------------------------------------------------------------------
1589      !! Function : w2 (see Gaussian function psi in Tsamados et al 2013)
1590      !!-------------------------------------------------------------------
1591      REAL(wp), INTENT(in   ) ::   pa
1592      REAL(wp) ::   w2
1593      !!-------------------------------------------------------------------
1594
1595      w2 = -    6670.68911883_wp &
1596       &   +   70222.33061536_wp*pa &
1597       &   -  314871.71525448_wp*pa*pa &
1598       &   +  779570.02793492_wp*pa*pa*pa &
1599       &   - 1151098.82436864_wp*pa*pa*pa*pa &
1600       &   + 1013896.59464498_wp*pa*pa*pa*pa*pa &
1601       &   -  493379.44906738_wp*pa*pa*pa*pa*pa*pa &
1602       &   +  102356.55151800_wp*pa*pa*pa*pa*pa*pa*pa
1603
1604   END FUNCTION w2
1605
1606   SUBROUTINE all_skr_sks( px, py, pz, allsk )
1607      REAL(wp), INTENT(in   ) ::   px,py,pz
1608      REAL(wp), INTENT(out  ) ::   allsk(6)
1609
1610      REAL(wp) ::   zs12r0, zs21r0
1611      REAL(wp) ::   zs12s0, zs21s0
1612
1613      REAL(wp) ::   zpih
1614      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
1615      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
1616      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22
1617      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22
1618      REAL(wp) ::   zd11, zd12, zd22
1619      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2
1620      REAL(wp) ::   zHen1t2, zHen2t1
1621      !!-------------------------------------------------------------------
1622
1623      zpih = 0.5_wp*rpi
1624
1625      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi)
1626      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi)
1627      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi)
1628      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi)
1629      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi)
1630      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi)
1631      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi)
1632      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi)
1633      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi)
1634      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi)
1635      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi)
1636      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi)
1637      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi)
1638      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi)
1639      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi)
1640      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi)
1641   ! In expression of tensor d, with this formulatin d(x)=-d(x+pi)
1642   ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else
1643   ! x=x-pi/2
1644      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))
1645      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))
1646      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))
1647      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22
1648      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22
1649      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22
1650
1651      IF (-zIIn1t2>=rsmall) THEN
1652      zHen1t2 = 1._wp
1653      ELSE
1654      zHen1t2 = 0._wp
1655      ENDIF
1656
1657      IF (-zIIn2t1>=rsmall) THEN
1658      zHen2t1 = 1._wp
1659      ELSE
1660      zHen2t1 = 0._wp
1661      ENDIF
1662
1663      !!-------------------------------------------------------------------
1664      !! Function : s11kr
1665      !!-------------------------------------------------------------------
1666      allsk(1) = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11)
1667      !!-------------------------------------------------------------------
1668      !! Function : s12kr
1669      !!-------------------------------------------------------------------
1670      zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12)
1671      zs21r0 = (- zHen1t2 * zn1t2i21 - zHen2t1 * zn2t1i21)
1672      allsk(2)=0.5_wp*(zs12r0+zs21r0)
1673      !!-------------------------------------------------------------------
1674      !! Function : s22kr
1675      !!-------------------------------------------------------------------
1676      allsk(3) = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22)
1677      !!-------------------------------------------------------------------
1678      !! Function : s11ks
1679      !!-------------------------------------------------------------------
1680
1681      allsk(4) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11)
1682      !!-------------------------------------------------------------------
1683      !! Function : s12ks
1684      !!-------------------------------------------------------------------
1685      zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12)
1686      zs21s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i21 + zHen2t1 * zt2t1i21)
1687      allsk(5)=0.5_wp*(zs12s0+zs21s0)
1688      !!-------------------------------------------------------------------
1689      !! Function : s22ks
1690      !!-------------------------------------------------------------------
1691      allsk(6) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22)
1692   END SUBROUTINE all_skr_sks
1693
1694#else
1695   !!----------------------------------------------------------------------
1696   !!   Default option         Empty module           NO SI3 sea-ice model
1697   !!----------------------------------------------------------------------
1698   USE par_kind
1699   USE lib_mpp
1700   CONTAINS
1701   SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12, prdg_conv )
1702      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
1703      INTEGER                 , INTENT(in   ) ::   Kmm                                   ! ocean time level index
1704      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pstress1_i, pstress2_i, pstress12_i   !
1705      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pshear_i  , pdivu_i   , pdelta_i      !
1706      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   paniso_11 , paniso_12                 ! structure tensor components
1707      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   prdg_conv                             ! for ridging
1708      CALL ctl_stop('EAP rheology is currently dsabled due to issues with AGRIF preprocessing')
1709   END SUBROUTINE ice_dyn_rhg_eap
1710   SUBROUTINE rhg_eap_rst( cdrw, kt )
1711      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
1712      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
1713   END SUBROUTINE rhg_eap_rst
1714#endif
1715
1716   !!==============================================================================
1717END MODULE icedyn_rhg_eap
Note: See TracBrowser for help on using the repository browser.