source: NEMO/trunk/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90 @ 14117

Last change on this file since 14117 was 14117, checked in by acc, 11 months ago

Re-working of icedyn_rhg_eap.F90 for better adhesion to coding standards and a reduction in the number of functions. AGRIF is no longer confused by this routine so the bypass introduced at revision 14020 has been removed. May need to force a clean remake of AGRIF-based SETTE tests to pick up changes. No changes in actual SETTE results.

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