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 @ 14433

Last change on this file since 14433 was 14433, checked in by smasson, 3 years ago

trunk: merge dev_r14312_MPI_Interface into the trunk, #2598

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