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

Last change on this file since 12261 was 12261, checked in by stefryn, 22 months ago

add yield surface output and include effect on ridging

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