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

Last change on this file since 13113 was 13113, checked in by stefryn, 3 months ago

added structure tensor evolution + some changes for NEMO coding rules, to be continued

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