source: NEMO/trunk/src/ICE/icedyn_rhg_evp.F90 @ 9660

Last change on this file since 9660 was 9660, checked in by clem, 3 years ago

remove boundary conditions for ice dynamics. It makes the ice rheology numerically unstable. However the effect of this removal on sea ice boundary conditions is still unclear

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