source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_rhg_evp.F90 @ 10402

Last change on this file since 10402 was 10402, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: no more need of lk_mpp for mpp_sum/max/min, see #2133

  • Property svn:keywords set to Id
File size: 55.2 KB
Line 
1MODULE icedyn_rhg_evp
2   !!======================================================================
3   !!                     ***  MODULE  icedyn_rhg_evp  ***
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   !!----------------------------------------------------------------------
16#if defined key_si3
17   !!----------------------------------------------------------------------
18   !!   'key_si3'                                       SI3 sea-ice model
19   !!----------------------------------------------------------------------
20   !!   ice_dyn_rhg_evp : computes ice velocities from EVP rheology
21   !!   rhg_evp_rst     : read/write EVP fields in ice restart
22   !!----------------------------------------------------------------------
23   USE phycst         ! Physical constant
24   USE dom_oce        ! Ocean domain
25   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
26   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
27   USE ice            ! sea-ice: ice variables
28   USE icevar         ! ice_var_sshdyn
29   USE icedyn_rdgrft  ! sea-ice: ice strength
30   USE bdy_oce , ONLY : ln_bdy 
31   USE bdyice 
32#if defined key_agrif
33   USE agrif_ice_interp
34#endif
35   !
36   USE in_out_manager ! I/O manager
37   USE iom            ! I/O manager library
38   USE lib_mpp        ! MPP library
39   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
40   USE lbclnk         ! lateral boundary conditions (or mpp links)
41   USE prtctl         ! Print control
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   ice_dyn_rhg_evp   ! called by icedyn_rhg.F90
47   PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
48
49   !! * Substitutions
50#  include "vectopt_loop_substitute.h90"
51   !!----------------------------------------------------------------------
52   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
53   !! $Id$
54   !! Software governed by the CeCILL license (see ./LICENSE)
55   !!----------------------------------------------------------------------
56CONTAINS
57
58   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
59      !!-------------------------------------------------------------------
60      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  ***
61      !!                             EVP-C-grid
62      !!
63      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
64      !!  stress and sea-surface slope. Ice-ice interaction is described by
65      !!  a non-linear elasto-viscous-plastic (EVP) law including shear
66      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).   
67      !!
68      !!  The points in the C-grid look like this, dear reader
69      !!
70      !!                              (ji,jj)
71      !!                                 |
72      !!                                 |
73      !!                      (ji-1,jj)  |  (ji,jj)
74      !!                             ---------   
75      !!                            |         |
76      !!                            | (ji,jj) |------(ji,jj)
77      !!                            |         |
78      !!                             ---------   
79      !!                     (ji-1,jj-1)     (ji,jj-1)
80      !!
81      !! ** Inputs  : - wind forcing (stress), oceanic currents
82      !!                ice total volume (vt_i) per unit area
83      !!                snow total volume (vt_s) per unit area
84      !!
85      !! ** Action  : - compute u_ice, v_ice : the components of the
86      !!                sea-ice velocity vector
87      !!              - compute delta_i, shear_i, divu_i, which are inputs
88      !!                of the ice thickness distribution
89      !!
90      !! ** Steps   : 0) compute mask at F point
91      !!              1) Compute ice snow mass, ice strength
92      !!              2) Compute wind, oceanic stresses, mass terms and
93      !!                 coriolis terms of the momentum equation
94      !!              3) Solve the momentum equation (iterative procedure)
95      !!              4) Recompute delta, shear and divergence
96      !!                 (which are inputs of the ITD) & store stress
97      !!                 for the next time step
98      !!              5) Diagnostics including charge ellipse
99      !!
100      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
101      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
102      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
103      !!              (i.e. more stable and better convergence)
104      !!
105      !! References : Hunke and Dukowicz, JPO97
106      !!              Bouillon et al., Ocean Modelling 2009
107      !!              Bouillon et al., Ocean Modelling 2013
108      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
109      !!-------------------------------------------------------------------
110      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
111      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
112      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
113      !!
114      INTEGER ::   ji, jj       ! dummy loop indices
115      INTEGER ::   jter         ! local integers
116      !
117      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio
118      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
119      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
120      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
121      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass
122      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
123      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                           ! temporary scalars
124      !
125      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
126      REAL(wp) ::   zintb, zintn                                        ! dummy argument
127      REAL(wp) ::   zfac_x, zfac_y
128      REAL(wp) ::   zshear, zdum1, zdum2
129      !
130      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
131      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
132      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
133      !
134      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points
135      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
136      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points
137      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
138      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
139      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
140      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
141      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
142      !
143      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
144      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
145!!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
146      REAL(wp), DIMENSION(jpi,jpj) ::   zssh_lead_m                     ! array used for the calculation of ice surface slope:
147      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
148      !                                                                 !    ice top surface if ice is embedded   
149      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
150      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
151      !
152      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
153      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
154      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
155
156      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
157      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
158      !! --- diags
159      REAL(wp), DIMENSION(jpi,jpj) ::   zswi
160      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
161      !! --- SIMIP diags
162      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
163      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice
164      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2)
165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2)
166      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2)
167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2)
168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2)
169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2)
170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress
171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress 
172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
174      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
175      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
176      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
177      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)     
178      !!-------------------------------------------------------------------
179
180      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology'
181      !
182!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
183      !------------------------------------------------------------------------------!
184      ! 0) mask at F points for the ice
185      !------------------------------------------------------------------------------!
186      ! ocean/land mask
187      DO jj = 1, jpjm1
188         DO ji = 1, jpim1      ! NO vector opt.
189            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
190         END DO
191      END DO
192      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
193
194      ! Lateral boundary conditions on velocity (modify zfmask)
195      zwf(:,:) = zfmask(:,:)
196      DO jj = 2, jpjm1
197         DO ji = fs_2, fs_jpim1   ! vector opt.
198            IF( zfmask(ji,jj) == 0._wp ) THEN
199               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
200            ENDIF
201         END DO
202      END DO
203      DO jj = 2, jpjm1
204         IF( zfmask(1,jj) == 0._wp ) THEN
205            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
206         ENDIF
207         IF( zfmask(jpi,jj) == 0._wp ) THEN
208            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
209         ENDIF
210      END DO
211      DO ji = 2, jpim1
212         IF( zfmask(ji,1) == 0._wp ) THEN
213            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
214         ENDIF
215         IF( zfmask(ji,jpj) == 0._wp ) THEN
216            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
217         ENDIF
218      END DO
219      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
220
221      !------------------------------------------------------------------------------!
222      ! 1) define some variables and initialize arrays
223      !------------------------------------------------------------------------------!
224      zrhoco = rau0 * rn_cio 
225
226      ! ecc2: square of yield ellipse eccenticrity
227      ecc2    = rn_ecc * rn_ecc
228      z1_ecc2 = 1._wp / ecc2
229
230      ! Time step for subcycling
231      zdtevp   = rdt_ice / REAL( nn_nevp )
232      z1_dtevp = 1._wp / zdtevp
233
234      ! alpha parameters (Bouillon 2009)
235      IF( .NOT. ln_aEVP ) THEN
236         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
237         zalph2 = zalph1 * z1_ecc2
238
239         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
240         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
241      ENDIF
242         
243      ! Initialise stress tensor
244      zs1 (:,:) = pstress1_i (:,:) 
245      zs2 (:,:) = pstress2_i (:,:)
246      zs12(:,:) = pstress12_i(:,:)
247
248      ! Ice strength
249      CALL ice_strength
250
251      ! scale factors
252      DO jj = 2, jpjm1
253         DO ji = fs_2, fs_jpim1
254            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
255            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
256         END DO
257      END DO
258           
259      !
260      !------------------------------------------------------------------------------!
261      ! 2) Wind / ocean stress, mass terms, coriolis terms
262      !------------------------------------------------------------------------------!
263
264      !== embedded sea ice: compute representative ice top surface      ==!
265      !== non-embedded sea ice: use ocean surface for slope calculation ==!
266      zssh_lead_m(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
267
268      DO jj = 2, jpjm1
269         DO ji = fs_2, fs_jpim1
270
271            ! ice fraction at U-V points
272            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)
273            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)
274
275            ! Ice/snow mass at U-V points
276            zm1 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
277            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
278            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * vt_i(ji  ,jj+1) )
279            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
280            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
281
282            ! Ocean currents at U-V points
283            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
284               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
285           
286            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
287               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
288
289            ! Coriolis at T points (m*f)
290            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
291
292            ! dt/m at T points (for alpha and beta coefficients)
293            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
294           
295            ! m/dt
296            zmU_t(ji,jj)    = zmassU * z1_dtevp
297            zmV_t(ji,jj)    = zmassV * z1_dtevp
298           
299            ! Drag ice-atm.
300            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
301            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
302
303            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
304            zspgU(ji,jj)    = - zmassU * grav * ( zssh_lead_m(ji+1,jj) - zssh_lead_m(ji,jj) ) * r1_e1u(ji,jj)
305            zspgV(ji,jj)    = - zmassV * grav * ( zssh_lead_m(ji,jj+1) - zssh_lead_m(ji,jj) ) * r1_e2v(ji,jj)
306
307            ! masks
308            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
309            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
310
311            ! switches
312            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
313            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
314
315         END DO
316      END DO
317      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zmf, 'T', 1., zdt_m, 'T', 1. )
318      !
319      !------------------------------------------------------------------------------!
320      ! 3) Solution of the momentum equation, iterative procedure
321      !------------------------------------------------------------------------------!
322      !
323      !                                               !----------------------!
324      DO jter = 1 , nn_nevp                           !    loop over jter    !
325         !                                            !----------------------!       
326         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
327         !
328!!$         IF(ln_ctl) THEN   ! Convergence test
329!!$            DO jj = 1, jpjm1
330!!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
331!!$               zv_ice(:,jj) = v_ice(:,jj)
332!!$            END DO
333!!$         ENDIF
334
335         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
336         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
337            DO ji = 1, jpim1
338
339               ! shear at F points
340               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)   &
341                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
342                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
343
344            END DO
345         END DO
346         CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1. )
347
348         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
349            DO ji = 2, jpi ! no vector loop
350
351               ! shear**2 at T points (doc eq. A16)
352               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
353                  &   + 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)  &
354                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
355             
356               ! divergence at T points
357               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
358                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
359                  &    ) * r1_e1e2t(ji,jj)
360               zdiv2 = zdiv * zdiv
361               
362               ! tension at T points
363               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)   &
364                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
365                  &   ) * r1_e1e2t(ji,jj)
366               zdt2 = zdt * zdt
367               
368               ! delta at T points
369               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
370
371               ! P/delta at T points
372               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
373
374               ! alpha & beta for aEVP
375               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
376               !   alpha = beta = sqrt(4*gamma)
377               IF( ln_aEVP ) THEN
378                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
379                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
380                  zalph2   = zalph1
381                  z1_alph2 = z1_alph1
382               ENDIF
383               
384               ! stress at T points
385               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
386               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
387             
388            END DO
389         END DO
390         CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1. )
391
392         DO jj = 1, jpjm1
393            DO ji = 1, jpim1
394
395               ! alpha & beta for aEVP
396               IF( ln_aEVP ) THEN
397                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
398                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
399                  zbeta(ji,jj) = zalph2
400               ENDIF
401               
402               ! P/delta at F points
403               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
404               
405               ! stress at F points
406               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
407
408            END DO
409         END DO
410
411         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
412         DO jj = 2, jpjm1
413            DO ji = fs_2, fs_jpim1               
414               !                   !--- U points
415               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
416                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
417                  &                    ) * r1_e2u(ji,jj)                                                                      &
418                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
419                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
420                  &                  ) * r1_e1e2u(ji,jj)
421               !
422               !                !--- V points
423               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
424                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
425                  &                    ) * r1_e1v(ji,jj)                                                                      &
426                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
427                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
428                  &                  ) * r1_e1e2v(ji,jj)
429               !
430               !                !--- u_ice at V point
431               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
432                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
433               !
434               !                !--- v_ice at U point
435               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
436                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
437            END DO
438         END DO
439         !
440         ! --- Computation of ice velocity --- !
441         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
442         !  Bouillon et al. 2009 (eq 34-35) => stable
443         IF( MOD(jter,2) == 0 ) THEN ! even iterations
444            !
445            DO jj = 2, jpjm1
446               DO ji = fs_2, fs_jpim1
447                  !                 !--- tau_io/(v_oce - v_ice)
448                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
449                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
450                  !                 !--- Ocean-to-Ice stress
451                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
452                  !
453                  !                 !--- tau_bottom/v_ice
454                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
455                  zTauB = - tau_icebfr(ji,jj) / zvel
456                  !
457                  !                 !--- Coriolis at V-points (energy conserving formulation)
458                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
459                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
460                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
461                  !
462                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
463                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
464                  !
465                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
466                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
467                  !
468                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
469                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
470                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
471                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
472                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
473                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
474                     &           ) * zmaskV(ji,jj)
475                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
476                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
477                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
478                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
479                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
480                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
481                     &            ) * zmaskV(ji,jj)
482                  ENDIF
483               END DO
484            END DO
485            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
486            !
487#if defined key_agrif
488!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
489            CALL agrif_interp_ice( 'V' )
490#endif
491            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
492            !
493            DO jj = 2, jpjm1
494               DO ji = fs_2, fs_jpim1         
495                  !                 !--- tau_io/(u_oce - u_ice)
496                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
497                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
498                  !                 !--- Ocean-to-Ice stress
499                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
500                  !
501                  !                 !--- tau_bottom/u_ice
502                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
503                  zTauB = - tau_icebfr(ji,jj) / zvel
504                  !
505                  !                 !--- Coriolis at U-points (energy conserving formulation)
506                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
507                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
508                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
509                  !
510                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
511                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
512                  !
513                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
514                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
515                  !
516                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
517                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
518                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
519                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
520                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
521                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
522                     &            ) * zmaskU(ji,jj)
523                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
524                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
525                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
526                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
527                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
528                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
529                     &            ) * zmaskU(ji,jj)
530                  ENDIF
531               END DO
532            END DO
533            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
534            !
535#if defined key_agrif
536!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
537            CALL agrif_interp_ice( 'U' )
538#endif
539            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
540            !
541         ELSE ! odd iterations
542            !
543            DO jj = 2, jpjm1
544               DO ji = fs_2, fs_jpim1
545                  !                 !--- tau_io/(u_oce - u_ice)
546                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
547                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
548                  !                 !--- Ocean-to-Ice stress
549                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
550                  !
551                  !                 !--- tau_bottom/u_ice
552                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
553                  zTauB = - tau_icebfr(ji,jj) / zvel
554                  !
555                  !                 !--- Coriolis at U-points (energy conserving formulation)
556                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
557                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
558                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
559                  !
560                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
561                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
562                  !
563                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
564                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
565                  !
566                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
567                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity
568                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
569                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
570                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0
571                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin
572                     &            ) * zmaskU(ji,jj)
573                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
574                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity
575                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
576                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
577                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
578                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
579                     &            ) * zmaskU(ji,jj)
580                  ENDIF
581               END DO
582            END DO
583            CALL lbc_lnk( 'icedyn_rhg_evp', u_ice, 'U', -1. )
584            !
585#if defined key_agrif
586!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
587            CALL agrif_interp_ice( 'U' )
588#endif
589            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
590            !
591            DO jj = 2, jpjm1
592               DO ji = fs_2, fs_jpim1
593                  !                 !--- tau_io/(v_oce - v_ice)
594                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
595                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
596                  !                 !--- Ocean-to-Ice stress
597                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
598                  !
599                  !                 !--- tau_bottom/v_ice
600                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
601                  ztauB = - tau_icebfr(ji,jj) / zvel
602                  !
603                  !                 !--- Coriolis at v-points (energy conserving formulation)
604                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
605                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
606                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
607                  !
608                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
609                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
610                  !
611                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction
612                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
613                  !
614                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
615                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity
616                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
617                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
618                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0
619                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin
620                     &           ) * zmaskV(ji,jj)
621                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
622                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity
623                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
624                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast
625                     &              + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
626                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin
627                     &            ) * zmaskV(ji,jj)
628                  ENDIF
629               END DO
630            END DO
631            CALL lbc_lnk( 'icedyn_rhg_evp', v_ice, 'V', -1. )
632            !
633#if defined key_agrif
634!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
635            CALL agrif_interp_ice( 'V' )
636#endif
637            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
638            !
639         ENDIF
640
641!!$         IF(ln_ctl) THEN   ! Convergence test
642!!$            DO jj = 2 , jpjm1
643!!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
644!!$            END DO
645!!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
646!!$            CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
647!!$         ENDIF
648         !
649         !                                                ! ==================== !
650      END DO                                              !  end loop over jter  !
651      !                                                   ! ==================== !
652      !
653      !------------------------------------------------------------------------------!
654      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
655      !------------------------------------------------------------------------------!
656      DO jj = 1, jpjm1
657         DO ji = 1, jpim1
658
659            ! shear at F points
660            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)   &
661               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
662               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
663
664         END DO
665      END DO           
666     
667      DO jj = 2, jpjm1
668         DO ji = 2, jpim1 ! no vector loop
669           
670            ! tension**2 at T points
671            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)   &
672               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
673               &   ) * r1_e1e2t(ji,jj)
674            zdt2 = zdt * zdt
675           
676            ! shear**2 at T points (doc eq. A16)
677            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
678               &   + 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)  &
679               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
680           
681            ! shear at T points
682            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
683
684            ! divergence at T points
685            pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
686               &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
687               &             ) * r1_e1e2t(ji,jj)
688           
689            ! delta at T points
690            zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 
691            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
692            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
693
694         END DO
695      END DO
696      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
697     
698      ! --- Store the stress tensor for the next time step --- !
699      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
700      pstress1_i (:,:) = zs1 (:,:)
701      pstress2_i (:,:) = zs2 (:,:)
702      pstress12_i(:,:) = zs12(:,:)
703      !
704
705      !------------------------------------------------------------------------------!
706      ! 5) diagnostics
707      !------------------------------------------------------------------------------!
708      DO jj = 1, jpj
709         DO ji = 1, jpi
710            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
711         END DO
712      END DO
713
714      ! --- divergence, shear and strength --- !
715      IF( iom_use('icediv') )   CALL iom_put( "icediv" , pdivu_i (:,:) * zswi(:,:) )   ! divergence
716      IF( iom_use('iceshe') )   CALL iom_put( "iceshe" , pshear_i(:,:) * zswi(:,:) )   ! shear
717      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
718
719      ! --- charge ellipse --- !
720      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
721         !
722         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
723         !         
724         DO jj = 2, jpjm1
725            DO ji = 2, jpim1
726               zdum1 = ( zswi(ji-1,jj) * pstress12_i(ji-1,jj) + zswi(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
727                  &      zswi(ji  ,jj) * pstress12_i(ji  ,jj) + zswi(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
728                  &    / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
729
730               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress 
731
732               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
733
734!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
735!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
736!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
737!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
738               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
739               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
740               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
741            END DO
742         END DO
743         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
744         !
745         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 )
746         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 )
747         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
748         !
749         DEALLOCATE( zsig1 , zsig2 , zsig3 )
750      ENDIF
751     
752      ! --- SIMIP --- !
753      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. &
754         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. &
755         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. &
756         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN
757
758         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  &
759            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  &
760            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  &
761            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) )
762         
763         DO jj = 2, jpjm1
764            DO ji = 2, jpim1
765               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
766               
767               ! Stress tensor invariants (normal and shear stress N/m)
768               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
769               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
770               
771               ! Stress terms of the momentum equation (N/m2)
772               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch     ! sea surface slope stress term
773               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
774               
775               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch     ! Coriolis stress term
776               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch
777               
778               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch     ! internal stress term
779               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
780               
781               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
782               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
783               
784               ! 2D ice mass, snow mass, area transport arrays (X, Y)
785               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
786               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
787               
788               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
789               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
790               
791               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
792               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
793               
794               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
795               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
796               
797            END DO
798         END DO
799         
800         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
801            &                zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
802            &                zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
803            &                zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
804                 
805         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zdiag_utau_oi  , 'U', -1., zdiag_vtau_oi  , 'V', -1.,   &
806            &                zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   &
807            &                zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   &
808            &                zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    )
809         
810         IF( iom_use('normstr' ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress
811         IF( iom_use('sheastr' ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress
812         IF( iom_use('dssh_dx' ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x)
813         IF( iom_use('dssh_dy' ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y)
814         IF( iom_use('corstrx' ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x)
815         IF( iom_use('corstry' ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y)
816         IF( iom_use('intstrx' ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x)
817         IF( iom_use('intstry' ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y)
818         IF( iom_use('utau_oi' ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x)
819         IF( iom_use('vtau_oi' ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y)
820         IF( iom_use('xmtrpice') )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s)
821         IF( iom_use('ymtrpice') )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport
822         IF( iom_use('xmtrpsnw') )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s)
823         IF( iom_use('ymtrpsnw') )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport
824         IF( iom_use('xatrp'   ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport
825         IF( iom_use('yatrp'   ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport
826
827         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  &
828            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  &
829            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  &
830            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     )
831
832      ENDIF
833      !
834   END SUBROUTINE ice_dyn_rhg_evp
835
836
837   SUBROUTINE rhg_evp_rst( cdrw, kt )
838      !!---------------------------------------------------------------------
839      !!                   ***  ROUTINE rhg_evp_rst  ***
840      !!                     
841      !! ** Purpose :   Read or write RHG file in restart file
842      !!
843      !! ** Method  :   use of IOM library
844      !!----------------------------------------------------------------------
845      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
846      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
847      !
848      INTEGER  ::   iter            ! local integer
849      INTEGER  ::   id1, id2, id3   ! local integers
850      !!----------------------------------------------------------------------
851      !
852      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
853         !                                   ! ---------------
854         IF( ln_rstart ) THEN                   !* Read the restart file
855            !
856            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
857            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
858            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
859            !
860            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
861               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
862               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
863               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
864            ELSE                                     ! start rheology from rest
865               IF(lwp) WRITE(numout,*)
866               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
867               stress1_i (:,:) = 0._wp
868               stress2_i (:,:) = 0._wp
869               stress12_i(:,:) = 0._wp
870            ENDIF
871         ELSE                                   !* Start from rest
872            IF(lwp) WRITE(numout,*)
873            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
874            stress1_i (:,:) = 0._wp
875            stress2_i (:,:) = 0._wp
876            stress12_i(:,:) = 0._wp
877         ENDIF
878         !
879      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
880         !                                   ! -------------------
881         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
882         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
883         !
884         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
885         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
886         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
887         !
888      ENDIF
889      !
890   END SUBROUTINE rhg_evp_rst
891
892#else
893   !!----------------------------------------------------------------------
894   !!   Default option         Empty module           NO SI3 sea-ice model
895   !!----------------------------------------------------------------------
896#endif
897
898   !!==============================================================================
899END MODULE icedyn_rhg_evp
Note: See TracBrowser for help on using the repository browser.