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

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

final variable updates

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