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

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

source: NEMO/branches/2020/tickets_icb_1900/src/ICE/icedyn_rhg_eap.F90 @ 14016

Last change on this file since 14016 was 14016, checked in by mathiot, 4 years ago

ticket 1900: upgrade to trunk@r14015 (head trunk at 16h27)

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