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

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

source: NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/ICE/icedyn_rhg_eap.F90 @ 14015

Last change on this file since 14015 was 14015, checked in by rlod, 4 years ago

upgrade branch to trunk@14010

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