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

Last change on this file since 13105 was 13105, checked in by stefryn, 4 months ago

small changes to calls

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